ABSTRACT 

BAJDICH, MICHAL. Generalized Pairing Wave Functions and Nodal Properties 

for Electronic Structure Quantum Monte Carlo. (Under the direction of Prof. Lubos Mitas.) 

The quantum Monte Carlo (QMC) is one of the most promising many-body elec- 
tronic structure approaches. It employs stochastic techniques for solving the stationary 
Schrodinger equation and for evaluation of expectation values. The key advantage of QMC 
is its capability to use the explicitly correlated wave functions, which allow the study of 
many-body effects beyond the reach of mean-field methods. The most important limit on 
QMC accuracy is the fixed-node approximation, which comes from necessity to circumvent 
the fermion sign problem. The size of resulting fixed-node errors depends on the quality of 
the nodes (the subset of position space where the wave function vanishes) of a used wave 
function. In this dissertation, we analyze the nodal properties of the existing fermionic wave 
functions and offer new types of variational wave functions with improved nodal structure. 

In the first part of this dissertation, we study the fermion nodes for spin-polarized 
states of a few-electron ions and molecules with s, p, d and / one-particle orbitals. We find 
exact nodes for some cases of two electron atomic and molecular states and also the first 
exact node for the three-electron atomic system in 4 5(p 3 ) state using appropriate coordinate 
maps and wave function symmetries. We analyze the cases of nodes for larger number of 
electrons in the Hartree-Fock approximation and for some cases we find transformations 
for projecting the high-dimensional nodal manifolds into 3D space. The nodal topologies 
and other properties are studied using these projections. Finally, for two specific cases of 
spin-unpolarized states, we show how correlations reduce the nodal structure to only two 
maximal nodal cells. 

In the second part, we investigate several types of trial wave functions with pairing 
orbitals and their nodal properties in the fixed-node quantum Monte Carlo. Using a set of 
first row atoms and molecules we find that the wave functions in the form of single Pfaffian 
provide very consistent and systematic behavior in recovering the correlation energies on 
the level of 95%. In order to get beyond this limit we explore the possibilities of expanding 
the wave function in linear combinations of Pfaffians. We observe that molecular systems 
require much larger expansions than atomic systems and that the linear combinations of 
a few Pfaffians lead to rather small gains in correlation energy. Further, we test the wave 
function based on fully-antisymmetrized product of independent pair orbitals. Despite its 
seemingly large variational potential, we do not observe significant gains in correlation 



energy. Finally, we combine these developments with the recently proposed inhomogeneous 
backflow transformations. 
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distorted due to spin correlations into CI nodal surface (right) with just two 
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4.1 Correlation energies obtained by QMC methods with the different trial wave 
functions: VMC and fixed-node DMC with HF nodes (HF) and STU Pfaf- 
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gains over HF nodes for BCS and STU Pfaffian wave functions. The statis- 
tical error bars are of the symbol sizes or smaller. Except for the Be atom 
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Pfaffian nodes (orange/medium gray); and the nodes of the CI wave function 
(yellow/light gray) in two different views (upper and lower rows). The CI 
nodal surface is very close to the exact one (see text). The HF node clearly 
divides the space into four nodal cells while Pfaffian and CI wave functions 
partitioning leads to the minimal number of two nodal cells. The changes in 

the nodal topology occur on the appreciable spatial scale of the order of 1 a.u. [69] 

5.1 Slater- Jastrow (SJ), Pfaffian- Jastrow (PF) and Cl-Jastrow (CI) wave func- 
tions with backflow (BF) correlations for carbon pseudopotential (He core) 
atom (left) and dimer (right) tested in VMC and DMC methods. Upper 
figure: Percentages of correlation energy versus a number of backflow pa- 
rameters. Notation: 2B for all electron-nucleus and electron-electron terms, 
3B for all electron-electron-nucleus terms only and 23B for all terms together. 
Lower figure: Variance of the local energy versus a number of terms [TH] 

5.2 Optimized electron-nucleus x{ r ) an d electron-electron u(r) two-body back- 
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Chapter 1 



Introduction 

The properties of quantum chemical and condensed matter systems of our every- 
day world are determined by the laws of quantum physics, which have been known since 
the 1930s. The distribution and motion of electrons surrounding the atomic nucleus are 
described by time-dependent Schrodinger equation^ but solving it for many electrons is ex- 
tremely hard. The difficulty comes from presence of electron-electron interaction term, 
which makes it impossible to separate this many-body problem into set of one-electron 
problems. In the past, the methods for solving the Schrodinger equation were based on 
replacing the difficult interaction term by some effective one, designed "to capture" the 
essential physics. The great success of these theories is the proof of the genius of these 
approximations. In this dissertation, we present a quantum Monte Carlo (QMC) method, 
which incorporates the electron-electron term directly and is able to solve the many-body 
Schrodinger equation almost exactly. Currently, it is the only method, which treats the full 
many-body problem and scales up to large systems. 

The term "quantum Monte Carlo" covers several different stochastic techniques 

adapted to determine either the ground state or finite-temperature equilibrium properties. 

Further, we restrict our discussion only to the ground state electronic properties. The 

simplest method explained is the variational Monte Carlo (VMC), which uses the stochastic 

integration for evaluation of expectation values for chosen trial wave function. Its accuracy 

1 Quite often it is necessary to solve the Dirac equation for core electrons and use relativistic (spin-orbit) 
corrections elsewhere. 
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sorely depends on the quality of a used trial wave function. This drawback is partially 
removed by diffusion Monte Carlo (DMC), which projects out the ground-state component 
of the starting trial wave function. The only inaccuracy of DMC comes from the fermion-sign 
problem, the inability to directly sample the fermionic wave function. It can be circumvent 
by the fixed-node approximation, which enforces the nodes (the subset of position space 
where the wave function vanishes) of a trial wave function on the projected wave function, 
but introduces a small fixed-node error. 

Fixed-node QMC simulations are typically more computationally demanding than 
traditional independent-particle quantum chemistry techniques, but have been very effective 
in providing high accuracy results for many real systems such as molecules, clusters and 
solids. Typically, for cohesive and binding energies, band gaps, and other energy differences 
the agreement with experiments is within 1-3% [HI [9]. The computational cost of QMC 
increases as the cube with system size, making calculations with hundreds of electrons 
tractable; the large clusters |10] and solids jllj up to 1000 electrons have already been 
studied. 

Today, the fundamental problem of accuracy improvement of QMC simulations lies 
in the elimination or at least in the control of the fixed-node errors. To achieve this, we focus 
in this thesis on studying the structure and properties of fermionic wave functions, as well 
as on finding better approximations to their nodes. This is quite challenging, because the 
fermion nodes are complicated high-dimensional manifolds determined by the many-body 
effects. Despite this difficulty, we were able to discover the exact nodes for several high- 
symmetry cases and to analyze the nodal structure of many spin-polarized and unpolarized 
systems (I2|H5|. 

In the second part of the thesis, we provide the partial answer to the search for 
better approximations to the nodes of fermionic wave functions. We propose a generalized 
pairing trial wave function in the Pfaffian functional form, which leads us to accurate and 
compact description of the nodes and results in overall improvement in the accuracy of 

qmc nans]. 

The last part of the thesis deals with yet another way how to improve the nodes 
of wave functions. It proposes the fermion coordinate transformation of a backflow type |16| 
[T71 [TBI [T9| l20l [2T| [22] , which was recently demonstrated to work also for inhomogeneous sys- 
tems [23, 24 . We generalize its application to Pfaffian pairing wave functions and perform 
the first tests of this approach. 



1.1. The Many-Body Schrodinger Equation 
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This thesis is organized as follows: 

• The remainder of this chapter gives an overview of the techniques used to study the 
quantum mechanics of many-body and chemical systems. 

• The Chapter [2] follows with the description of the methodology behind the VMC and 
DMC methods. 

• Our study of the exact and approximate nodes of fermionic wave functions is summa- 
rized in the Ch. [3] 

• Chapter [4] discusses our calculations with generalized pairing wave functions. 

• Some preliminary results on the backflow corrected trial wave functions are the content 
of the Ch. 1 

• Finally, the last chapter concludes this dissertation. 

1.1 The Many-Body Schrodinger Equation 

One of the main challenges of condensed-matter physics and quantum chemistry 
is the accurate solution of the Schrodinger equation. Since the nuclei are about thousand 
times heavier than the electrons, the most common approach is to decouple the electronic 
and ionic degrees of freedom that is known as the Born- Oppenheimer approximation. The 
electronic part of non-relativistic Born-Oppenheimer Hamiltonian for quantum system of 
N e electrons in the presence of Nj nuclei in Hartree atomic units (h = m e = e = ^tt/cq = 1) 
is then given by 



where i and j indexes are summing over all electrons and / index is summing over all nuclei. 
A spectrum of states {& n } which diagonalizes the stationary Schrodinger equation 





n 




with Hamiltonian (1.1 ) is the solution of above many-body problem. 



1.2. Hartree-Fock Theory 
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The history of ab initio method:^] for electronic structure started soon after the 
invention of quantum mechanics by pioneering work of Heitler and London [26] , who calcu- 
lated the binding energy of H2. The first quantitative calculations of multi-electron systems 
where accomplished by Hartree (27] , and Hylleraas |28| 2!) . Fock (30] was initial in using the 
properly antisymmetrized wave function, the Slater (31] determinant, which established the 
Hartree-Fock (HF) theory. The HF theory replaces the hard problem of many interacting 
electrons with system of independent electrons in self-consistent field (SCF). 

In electronic structure of solids, the discoveries of the effective free electron theory 
together with Pauli exclusion principle [32 and the band theory of Bloch [33] were the 
first critical steps toward understanding crystals. In 1930s, the foundations for the basic 
classification of solids into metals, semiconductors and insulators were laid. Soon after, 
Wigner and Seitz [33] [35] performed the first quantitative calculation of electronic states of 
sodium metal. 

Today, the density functional theory (DFT) invented by Hohenberg and Kohn |36j 
and applied by Kohn and Sham [37] is the principal method for calculations of solids. 
Together with HF and post HF methods, they are very relevant to our discussion of quantum 
Monte Carlo (QMC), which relies on these methods mainly for its input. What follows is 
their brief overview. 



1.2 Hartree-Fock Theory 

Due to the Pauli exclusion principle, any solution to the stationary Schrddinger 
equation with Hamiltonian (1.1) has to be antisymmetric under exchange of any two elec- 
trons with the same spin as 

*(... .. ..) = -*(•• (1.3) 



For more details of the history of electronic structure, refer to a book by R. M. Martin 



1.3. Post Hartree-Fock Methods 
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The Hartree-Fock theory [301 131] uses the simplest antisymmetric wave function, the Slater 
determinant, 



VP(l,2,...,iV) 



^^(-1)^(1)^(2)...^^) 
iV. p 



(1.4) 



(p x (l) 0i(2) 

&(1) &(2) 

(^at(I) <^tv(2) 



£ 2 (7V) 
^w(iV) 



det^i(l),...,^(iV)], 



to approximate the state of N electron system. The (f>i(j) are one particle spin-orbitals, 
each of which is a product of spatial tp%(J) and spin afaj) orbitals. (Note that <pf(j) is 
independent of spin a in closed-shell cases. In open-shell systems, this assumption corre- 
sponds to spin-restricted Hartree-Fock approximation). If the Hamiltonian is independent of 



spin or is diagonal in spin basis a = {| j), | J.)}, the expectation value of Hamiltonian (1.1 ) 



with the wave function (1.4) is given by 
Ehf = 



<*(r) 



exi 



+ 



1 



E 



E 



<(r)dr 



dr dr' 



C(r)^(r')<(r)<(r') 



dr dr'. 



(1.5) 



Above expression for HF energy contains three different terms, the first being just sum 
of independent one-particle energies in the external potential of nuclei V ex t, the second is 
the direct contribution to Coulomb interaction also called the Hartree term and the last 
one is the exchange term (or Fock term). Note, that in the case of i = j (self-interaction 
contribution) the last two terms explicitly cancel each other. 



1.3 Post Hartree-Fock Methods 

The single-determinant HF theory contains proper antisymmetry, which introduces 
the effects related to the exchange. However, the full electron-electron Coulomb repulsion is 
only approximated by Hartree term. What is left out is referred to as electronic correlation. 
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The correlation energy is then defined as a difference of energies of exact wave function and 
the best HF wave function of the same state. Typically, the correlation energy constitutes 
only a fraction of total energy, but it accounts for large portion of cohesion and excitation 
energies. Obtaining the missing correlation is therefore the principal challenge of all modern 
electronic structure methods. The high accuracy is needed to access such problems as 
description of magnetism or superconductivity. 

The missing correlation in HF wave function can be accounted for by additional 
Slater determinants. There exist several methods, which can produce these multi-determinantal 
wave functions. The most frequently used are configuration interaction (CI), multi-configurational 
self-consistent field (MCSCF), and coupled cluster (CC) methods. The idea behind CI is 
to diagonalize the iV-electron Hamiltonian on the space of orthogonal Slater determinants. 
These determinants are typically constructed as multi-particle excitations from reference 
determinant (usually HF solution). If all the determinants are included, the full CI is in 
principle exact. However, the number of terms grows exponentially with N and therefore, 
we have to do limited expansions in practice. The most used is the configuration interac- 
tion with all the single and double excitations (CISD). The disadvantage of this approach 
is that the expansion converges very slowly in correlation energy (as y/~N) and creates the 
size- consistency problem (i.e., wrong scaling of total energy with N). 

The MCSCF is in some way a modification of truncated CI expansion, when the 
orbitals used for construction of determinants are also optimized together with determinan- 
tal weights. The optimization of all parameters is a difficult task and limits the number of 
determinants in the expansion. 

The size-consistency problem can be overcome by coupled cluster expansion [38 . 
All the excitations from the reference determinant are in principle included, but the coef- 
ficients of expansion are approximated and method is non-variational. The computational 
cost for CCSD (with singles and doubles) scales as N G and therefore constitutes a formidable 
challenge for larger systems. The commonly used CCSD(T) (CC with singles, doubles and 
approximated triples) produces the most accurate energies for small and intermediate sys- 
tems available at the present time, which in many cases serve as benchmarks for all other 
methods including QMC. For excellent treatment and overview of above quantum chemical 
methods we refer the reader to book by A. Szabo and N. S. Ostlund [39J. 
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1.4 Density Functional Theory 

Previous methods were examples of the wave function based theories. For density 
functional theory the primary object is the one-particle electron density. It is formally 
an exact method based on Hohenberg and Kohn theorem [36] that states that the ground 
state properties of a many-electron system can be obtained by minimizing the total energy 
functional 

E[p(r)] = J V ext (r)p(r) dr + F[p(r)}, (1.6) 

where the F[p(r)] is some unknown, universal functional of electronic density p(r). The total 
energy E has a minimum when the electronic density p(r) is equal to the exact electronic 
density in some external potential V ex t(r). 

Kohn and Sham [37] wrote the ansatz for the density in the terms of one-electron 
orbitals of an auxiliary non-interacting system as 

N 

p(r) = J>,(r)| 2 . (1.7) 

i 

The total energy of electronic system can be then expressed as 

#)] = -^E / ^(r)VV*(r)dr + J p(r)V ext (r) dr + \ jj drdr' + E xc [p{v)] , 

(1.8) 

where the first 2 terms represent the familiar energy of non- interacting system in the external 
potential, the third term is just the Hartree term and the rest is an unknown universal 
exchange- correlation functional of density. If the E xc [p(r)] would be precisely known, we 
would have one to one mapping between difficult many-electron system and this one-particle 
problem. 

The simplest approximation to exchange-correlation functional is the local density 
approximation (LDA), 

E xc [p(r)] = j ^(p(r))p(r)dr, (1.9) 

where the ext 9 is the exchange-correlation energy per electron in homogeneous electron gas. 
It is interesting to know that for practical application of the LDA, the correlation portion 
of this energy was taken from high-precision QMC calculation [ID] . Even today, the LDA 
is widely used in solids. 

In cases where the LDA is not accurate enough, it seem natural to express the 
e xc(p( r )> Vp(r)) as a function of the local density and its gradient. This was the essential 
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idea behind the generalized gradient approximation (GGA) |41l |4"2] 143] . which increased 
precision enabled wide-spread use in quantum chemical systems. 

There exist many flavors of density functional theory, e.g., hybrid-functionals |41| 
H2J SI], time-dependent DFT (TD-DFT) |45l S6J, extensions to non-local density func- 
tionals (averaged density approximation (ADA) [47 and weighted density approximation 
(WDA) |47j). or orbital-dependent functionals (e.g. self-interaction corrections (SIC) |1H] 
and LDA+U 09]). 

The best known density functional results are typically an order of magnitude less 
accurate than good QMC results. However, the computational cost is much more favor- 
able for DFT and consequently DFT is much more widely applied to variety of interesting 
applications in many fields of science and technology. 
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Chapter 2 

Quantum Monte Carlo Methods 

2.1 Introduction 

There are several Quantum Monte Carlo methods. They all apply a stochastic 
approach to find a solutions of a stationary Schrodinger equation of quantum systems. 
However, in this dissertation we will be restricted to description of only the variational 
and diffusion Monte Carlo methods. In the variational Monte Carlo the expectation values 
are evaluated via stochastic integration over 3N dimensional space. A variational theorem 
ensures that the expectation value of Hamiltonian with respect to given trial wave function 
is a true upper bound to the exact ground-state energy and sorely depends on the accuracy 
of the trial wave function. The second method, the diffusion Monte Carlo, removes some 
deficiency in the accuracy of trial wave function by employing imaginary-time projection 
of trial wave function onto the ground state. The nodal structure of trial wave function is 
enforced on the projected wave function by fixed node approximation. Hence, if the nodes 
of trial wave function were identical to the nodes of an exact one, we would know the ground 
state energy of a quantum system in polynomial time. 

2.1.1 Metropolis Sampling 

Most of quantum Monte Carlo calculations work with the Metropolis rejection 
algorithm [50] , in which a Markov process is constructed to generate a random walk through 
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the state space. (In VMC and DMC methods we walk in 3iV dimensional spatial space {R} 
and one fixed point R is also commonly referred to as walker.) The algorithm enables 
us to sample desired multi-dimensional probability distribution P(R) without any prior 
knowledge of its normalization. Given some transition rule P(R — > R') from R to R', 
which satisfies ergodicity and detailed balance, 

P(R-»R')P(R) = P(R' R)P(R'), (2.1) 

the desired probability density P(R) will converge to an equilibrium state given by 

P(R^R')P(R)dR = P(R). (2.2) 



/ 



The transition rule P(R — > R') can be further written as a product of the sampling distri- 
bution T(R — > R') (typically Gaussian distribution centered around R) and the probability 
of an acceptance ^4(R — > R') of the proposed step 

P(R-»R') = T(R->R')A(R-»R'). (2.3) 

The conditions of detailed balance are satisfied by choosing the A(R — > R') to be 

T(R -► R)P(R) 



A(R -> R') = min 



1. 



(2.4) 



T(R -► R')P(R) 

A Monte Carlo (MC) simulation initialized at some random state evolves toward an equi- 
librium. After the stabilization period we reach the point, where we can start to collect the 
statistics. The usual estimators for the mean, variance and error bars of the mean of some 
operator O(R) are given as 

M 

M 



m 

M 

( ^ )M = M^T^ (0(Rm) " (0)M)2 ' 



M _ a O 



>M 

where M represents a number of sampling points. Furthermore, the well known central 
limit theorem then ensures that 



hm (0) M = (O) (2.6) 

M^OO 



with statistical error bars going to zero as 1/y/M. 
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2.2 Variational Monte Carlo 

Variational Monte Carlo (VMC) performs stochastic evaluation of an expectation 
value of Hamiltonian H over trial wave function 'I'x'(R), which is a reasonably good ap- 
proximation of the true ground state. The is subject to several necessary conditions, 



which are addressed in Sec. 2.4 of this chapter. 

As follows from variational theorem, this expectation value provides rigorous upper 
bound on ground state energy Eq, 

/^(R)7flMR)dR J|*T(R)| 2 £L(R)dR 
EVMC ~ /tt$,(R)tt T (R)dR " /|* r (R)|*dR " °' ( ?) 

where in the second step we have introduced the term called local energy 

Wr(R) 

El(r) " ^(rT (2 - 8) 

evaluated over probability density (also called the importance function) 

v(Tt) _ l^(R)| 2 , . 

vm ~ ;i^(R)l 2 dR- (2 - 9) 

The estimator of Evmc is then given as 

1 - 



where M configurations are distributed according 'P(R) (2.9) via Metropolis algorithm. 

Trial wave function typically depends on a set of variational parameters, which 
can be optimized in order to achieve the minimum of Ey mc or the minimum of variance of 
local energy 

2 _ /|^r(R)| 2 (^L(R)-^Afc) 2 dR 
avMC ~ /|^ T (R)|2dR ' (2 - H) 

The variance o~ VMC is especially good function to minimize, since it is always positive and 
bounded from bellow (o~ VMC — ► as — * 3>o)- I n practice, there are other possible com- 
binations of EvmCi a vMC an d some other functions of which serve as good minimizers. 



We devote the Sec. 2.6 of this chapter to optimization of variational trial wave functions. 
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2.2.1 Correlated Sampling 



The technique of a correlated sampling exploits the fact that the variance of differ- 
ence of some correlated random variables X and Y decreases as their correlation increases. 
In QMC, the standard application of correlated sampling is to the VMC energy differences 
of two close (i.e., correlated) trial-wave functions, Vl/ip and . The energy difference can 



be written as 



E\ — E2 



J \*£\B)\*E£>(B.) dR J |^(R)| 2 ^(R) dR 

/l4 2) (R)p d R 

a)/ 



(2), 



,(2), 



J|*«(R)PdR 



4 1} (R) 



M 1} (R) 



J|*^(R)| 2 dR Ju>(R)|^(R)| 2 dR 



dR, 



(2.12) 



where in the second line we used re- weighting factor w(R) 



l4 2) (R)l 2 

|*W(R)I2' 



Therefore, we 



can obtain the above energy difference as a sum over weighted local energy differences. 
Correlated sampling in DMC, however, requires further modification to the Green's function 
(best known is the Filippi and Umrigar's approximation [51J. 



2.3 Diffusion Monte Carlo 

2.3.1 Imaginary Time Schrodinger Equation 

Diffusion Monte Carlo method belongs to a larger class of projection and Green's 
function MC methods. It exploits the imaginary time many-body Schrodinger equation 

d^fR t) 

} = -(H - £ T )*(R, r), (2.13) 
where ^(R, r) is the projected wave function with real variable r measuring the progress 



in imaginary time and Et is an energy offset. The equation (2.13) is formally similar to a 
diffusion equation. Its effect is to converge the initial wave function to the ground state. 
This can be easily seen by expanding ^(R, r) in the eigenstates {& n } of TL and projecting 
time r — > 00 

lim *(R,r) = lim V($ n |^(0)> e - ^"-^ $ n (R) (2.14) 

n 

= lim ($ |*(0)> e^^-^^oCR)- 
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Therefore, as far as there is non-zero overlap of the initial wave function with a ground state, 
we can always project out the later one. The parameter Et is adjusted during simulation to 
keep the amplitude of Vl/ constant. The initial wave function is typically analytically known 
[e.g. \P(R, 0) = ^t^R)] and the final projected wave function is represented via ensemble 
of walkers. Therefore, in this case the ^ not \l/ 2 is considered to be a population density of 
walkers. 

The diffusion algorithm is performed via Green's function 

G(R->R',t) = (Rje-^-^lR'), (2.15) 

which determines the wave function after projection time r as 

*(R,t) = yG(R^R , ,r)^(R / ,0)dR / . (2.16) 

The Green's function's short-time approximation |52j can be written using Trotter-Suzuki 
formula as 

G(R -> R',r) = (R|e- r ( f+ ^-^)|R') (2.17) 
_ e -r(V(R)-B)/2/ R , e -rf | R /\ e -T(V(R')-Br)/2 

« (2vrr)~ 3iV / 2 e X p 

where in the last line we used the solution for kinetic term as a Gaussian expanding in 
time. We can interpret the action of the short-time Green's function as a diffusion process 
(kinetic energy term) and branching or re-weighting process (potential-energy term) on the 
population of walkers. 



(R — R ) 2 

2r 



V(R) + V(R) 

exp -r ( E T 



2.3.2 Importance Sampling 

Recent DMC methods treat diffusion problem by introducing the importance sam- 
pling [5H (531 G3] with the mixed distribution /(R,r) = *r(R)*(R,r), where \I>t(R) is 
some trial wave function. The modified diffusion equation for above mixed distribution / 
is then 

d/( ^; ,r) = - ^V 2 /(R, t) + V • [v D (R)/(R, r)] (2.18) 
+ (E L (R)-E T )f(R,r), 
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where v_d(R) = Vln |^j-(R)| represents a drift velocity term and -Ex(R) is the local energy 



of ^j'(R) from Eq. (2.8). The short-time approximation to the Green's function is then 
correspondingly 

(R - R' - tv£)(R')) 2 " 



G{R — > R',r) ^(2ttt) 



-37V/2 



x exp 



6XP 2r 
/ E L (R) + E L (B!) 



(2.19) 



Ex 



V 2 
= Gd(R — >■ R', r) x G b (R->R',t). 

Here, we can again interpret the action of the short-time Green's function as a diffusion 
process Go (kinetic energy term) and branching/re- weighting process Gb controlled by 
difference of averaged local energies at R and R' to Et- This transformation has several 
consequences. First, the density of walkers is enhanced in the regions, where \&r(R) is 
large due to drift velocity v/)(R). Second, the re- weighting term contains the local energy 
£x(R) in the exponent, which for a good is close to a constant and therefore much 
better behaved than unbounded potential ^(R). 

In the actual calculation, we start with walkers distributed according to /(R, 0) = 
|^t(R)| 2 with the weights {w} of all walkers set to one. Then we perform MC step sampled 



from the kinetic part Go of Eq. (2.19). The step is accepted with Metropolis probability 

t G d (R-R)|*t(R')I 2 \ 



ACR -> R) = min 1, 

^'G (R-Rf r (R)P 

Consequently, we assign to each walker a new weight according to 



(2.20) 



W; 



G B (Ri ->Rj,T)u> f . 



(2.21) 



The advantage of this approach is that the DMC algorithm has essentially the VMC dynam- 



ics with a small time step with additional re-weighting [Eq. (2.21)]. Over the simulation, 
some weights will start to dominate over others. It is therefore necessary to control their 
population. One way is to take two walkers, one with a large weight w\ and another with 



a small weight W2- The first walker is then branched with probability 



Wl+W 2 



to two, both 



having an average weight of Wl + W2 . Hence the walker with a small weight gets killed. Af- 
ter an equilibration period, we can start to collect statistics needed for the calculation of 
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projected ground state energy 



r^*(R,r)W^ T (R)dR 

g °"° = £g. V(iU)«Mk)dR (2 ' 22) 
C/(R,T)Ei(R)dR 

= fun 7- — ; ; — ~ 

™ //(R,r)dR 

_ ZiWiE L (Ri) 
Ei ™i 

2.3.3 Fixed-Node Approximation 

For the bosonic systems, the ground state wave function (also ^(R, t) for that 
matter) has no nodes and can be considered to be positive everywhere. However, the 
fermionic wave function is naturally antisymmetric, and therefore its amplitude can no 
longer represent the sampling distribution for the walkers, which has to be positive every- 
where. We could try to overcome this problem in naive way by decomposing \&(R, r) into 
positive functions \£ + (R, r) and Vl/~(R, r) such that 

*(R, r) = ^ + (R, t) - ~(R, r) (2.23) 

However, since our diffusion equation is linear in VI/(R,t), both ^ + (R, r) and ^~(R, r) 
will converge to the same bosonic distribution, hence making the averaging of fermionic 



observables difficult (See Fig. 2.1). This is a well-known fermion sign problem. Up to day, 
no satisfactory solution has been found, however it can be circumvented by the fixed-node 
approximation |52[ 1551 156) [57] , which restricts the evolution of positive and negative walkers 
into separate regions defined by sign of trial wave function Vl/y(R). This is achieved by 
introducing a mixed distribution /(R, r) = ^ r r(R)^ r (R, t) in the importance sampling, so 
that we can extend the stochastic sampling of the walkers for fermionic wave function to 
the regions where 

/(R,r)>0 (2.24) 
and prohibiting it everywhere else. In other words, by restricting the Vl/(R, r) to have the 



same sign as \I/t(R) for all r we can readily satisfy Eq. (2.24). Therefore, the fixed-node 



approximation reduces the nonlocal antisymmetric condition for fermionic wave function to 



local boundary condition in Eq. (2.24). Furthermore, the drift-velocity term diverges at a 



node of ^ / t(R) = and acts like a repulsive potential, which stabilizes the simulation. The 
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x=0 




Figure 2.1: Imaginary time behavior of walker distributions \I /+ (red) and (blue). Upper 
figure: at the beginning of the diffusion (r = 0). Lower figure: After large enough imaginary- 
time evolution. The fermionic signal decays as — oc exp[— t(Eq — Eq)], where 
Eq — Eq is a difference of energies of fermionic to bosonic ground states. 

fixed-node DMC energy is still a rigorous upper bound to the true ground state energy |52| 
[57] . The resulting fixed-node error is proportional to the square of the nodal displacement 
error. 

2.3.4 Beyond Fixed-Node Approximation 

The fixed-node approximation was historically the first and the simplest way how 
to overcome the spurious fermion sign problem. Besides improving the nodes with better 
trial wave functions (see Chs. [4]and[5]) or attempts to understand them (see Ch. [3]), there 
have been several efforts to go beyond the fixed-node approximation. One of the first were 
ideas based on the cancellation of paired walkers with opposite signs [SS] • This method was 
subsequently improved to work with fully interacting ensemble of walkers [59] . Although the 
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method suffered from unfavorable scaling caused by computation of inter-walker distances, 
it was applied (among others) to Li atom ground state with good agreement with estimated 
exact energy. 

Ceperley and Alder [60 on the other hand developed the released-node method. 
Their method starts from an initial fixed-node distribution, but walkers are allowed to 
cross the nodes of a trial wave function. Each time a walker crosses the node it caries 
the additional weight in the form of db sign, which contributes to the released-node energy. 
This algorithm is in principle transient (the decay of fermion first excited state has to be 
faster then the decay of fermion ground state to bosonic state), the authors were able to 
successfully calculate the total energies of a number of small atoms and molecules. 

Subsequently, Anderson and Traynor [61] combined the best features of fixed-node, 
released-node and cancellation methods to algorithm, which employed improvements as 
relocation after node crossing, self-cancellations and multiple cancellations, and maximum 
use of symmetry in cancellations together with rigorous energy evaluation using importance 
sampling with trial wave functions. Among the most interesting systems, the authors 
calculated the excited state of H2 and the barrier height for the reaction H + H2 — > 
H2 + H [62]. However, the main obstacle to the application of this method to large systems 
is the computational cost associated with efficient annihilation of walkers. 

Among other alternative methods, let us also mention the method of adaptive nodes 
(see, e.g., Ref. |63j). The principle is based on representing the approximate nodal function 
directly from ensemble of walkers (e.g. Gaussian centered on each walker). This results 
in the adaptive description of the nodes that does not depend upon a priori knowledge 
of the wave function. The main difficulty of the method however comes from scaling at 
higher dimensions N, where the antisymmetry condition results in the calculation of AM 
permutations. 

The critical component of any method, which would solve the now famous fermion 
sign problem, is the polynomial scaling with the number of particles. There might exist 
approaches, which scale polynomially for some cases [63], but in general will scale expo- 
nentially. Indeed, recently has been shown [65] that the problem of quantum Monte Carlo 
simulation of a random coupling Ising model with fermion sign problem belongs into the 
class of nondeterministic polynomial (NP) hard problems. It is widely excepted that no NP 
problem can be solved in the polynomial (P) time, hence NP 7^ P. This finding thus implies 
that the general solution to the fermion sign problem in polynomial time is not possible. 
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2.4 Variational Trial Wave Functions 

The great advantage of QMC methods is the capability to use explicitly correlated 
trial wave functions. Their quality is the most important factor on final accuracy and 
efficiency of QMC. On the other hand, the repeated evaluation of (and V^r, ^7 2 ^t) 
is the most costly part of QMC calculation. Therefore, it is desirable to seek both highly 
accurate and efficient representations of 



2.4.1 Basic Properties 

Any reasonably accurate trial wave function for QMC has to obey several basic 
properties. Since is the approximate solution to the bound electronic system, we de- 
mand the / ^x^T, J ^>* T Ti^T and J ^* t H 2 ^t to exist. This is obvious, since the VMC 
expectation value and its variance would not be properly defined otherwise. In addition, if 
Hamiltonian does not explicitly depend on the magnetic field (has time-reversal symmetry) , 
$>T can be made real. 

Further, and VvE't has to be continuous, wherever the potential is finite. In 
addition, the continuity of the potential implies also the continuity of V 2 ^t- Important 
corollary is as approaches an exact solution, local energy becomes constant everywhere. 
Therefore, any singularity from the potential has to be properly canceled by an opposite 
singularity in the kinetic energy, which gives rise to Kato cusp conditions |66| [BT] . As the 
distance of electron to nucleus goes rj $ — ► 0, the potential energy divergence gets properly 
canceled when 



i dv T 

*T dm 



= -Zj. (2.25) 

r H =0 



The electron-nucleus cusp condition [Eq. (2.25)] is typically satisfied by the choice of proper 



orbitals or removed by the use of pseudopotentials (see Sec. 2.5). Similarly, as the electron- 
electron distance goes nj — ► 0, there is a divergence in electron-electron Coulomb potential. 
Let us for now assume that i and j electrons bear the same spin. If we write trial wave 
function as a product of an antisymmetric function ^A{fij) and a symmetric function 
^s{ r i,j) with respect to rij, then the divergent terms cancel if 



1 d^ s 



\- (2-26) 
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For two electrons with unlike spins, the situation is different, since the antisymmetric part 
^A(fij = 0) is generally non-zero. The unlike spins cusp condition therefore modifies to 



1 d^ s 



\ (2-27) 



*S d rij 

Finally, has to be antisymmetric with respect to exchange of both spin and spatial 
coordinates 

^ T (^R, P£) = (-1) P ^t(R, S), (2.28) 

where S = (a±, ...,<jn) are discrete spin variables with values ±i and P is an arbitrary 
permutation with sign equal to (— l) p . However, if the Hamiltonian does not contain spin- 
dependent terms, we will impose antisymmetry only with respect to interchanges between 
the electrons of the same spin. That is, the spatial antisymmetry is only required when 

PY. = S. (2.29) 

Then the spins of all electrons are being fixed with total spin projection equal to h(N^-N^). 
We can therefore label the first iVj particles as spin-up and the rest N — N\ as spin-down. 



Throughout this dissertation we will always assume that Eq. (2.29) holds and that the spin 
variable is factored out by the above labeling scheme. 



2.4.2 Trial Wave Functions Forms 

Almost all trial wave functions used in the QMC today are expressed as product of 
some antisymmetric part ^^(X) and symmetric exponential of Jastrow correlation factor 

tf T (R) = *x(X) x exp[£/ corr (R)], (2.30) 

where X = (xi, . . . ,xjv) represents some general quasi-particle coordinates dependent on 
all N electron positions R. For simplicity, let us assume X = R. We will come back to a 
more general case when we will be discussing the backflow transformation. 

Bellow, we will describe the particular forms for ^^(R), f7 corr (R) and X(R) as im- 
plemented in our QMC code QWALK [3 . Jastrow part [191 EZ1 EH1 [69] , explicitly dependent 
on electron-nucleus and electron-electron distances, 

Ucorr({nj},{ r ii},{ r ji}) = ^2x(ni) + ^u(r ii ) + ^ w(nj, r iI} r 5I ) (2.31) 

il i^j 
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is written as the sum of one-body (electron-nucleus), two-body (electron-electron) and three- 
body (electron-electron-nucleus) terms. The function of these terms is twofold. First is to 



describe the electron-electron cusp conditions [Eqs. (2.26) and (2.27)]. These are satisfied 
if 



-7 (Jj = <7,-, . 

4 3 (2.32) 

2 a i = ~ a ji 



and all the remaining functions have zero derivative. Second purpose is to incorporate 
the correlation effects not present in ^^(R). Expanded in the basis of one dimensional 
functions the components of the Jastrow factor take form 

x(r) = Y,4 n Mr), (2.33) 

k 

u(r) = Y,ci e b k (r), (2.34) 

k 

w{rij,rii,rji) = ^ c|f™ [a k (ru)ai(rji) + a fc (r j /)a i (r i /)]6 m (r ii ). (2.35) 

klm 

The above Jastrow factor has proved to be very efficient in describing correlation effects 
with minimal number of variational parameters. More details about the implemented basis 
functions can be found in App. |A]or in the QWalk documentation [3]. 

To propose a general fully antisymmetric wave function is a formidable chal- 
lenge given the scarcity of antisymmetric algebraic forms with known evaluation schemes. 
This dissertation deals with two of them, determinant and Pfaffian. The great success of 
quantum-chemical theories is based on the wave function constructed from Slater determi- 
nants written in spin-factorized form as 

*f ater (R) = J2m d}(1, . . . , JV T ) D\{N } + 1, . . . , N). (2.36) 

i 

Each spin-up and spin-down Slater determinant 

D](l, . . . , JV T ) = det[v?i T (l), . . . , <4 T (7V T )], (2.37) 
Df(JV T + 1, . . . , N) = det#J(JV T + 1), . . . , y\ (AO] (2.38) 

is a function of one-particle orbitals taken from (post-)HF or DFT methods. In the introduc- 
tory chapter we have touched upon some of the properties of \^j^ ater . The most important 
one is that the Slater determinants can form a complete set of antisymmetric functions and 
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are therefore in principle capable of describing the exact wave function. In practice, how- 
ever, we are limited to finite expansions. The missing correlation is then added by Jastrow 



is most commonly used type of trial wave function in QMC. Although very successful, the 
^T~ J has its limitations summarized bellow: 

(a) The large Cl-like expansions needed for high accuracy results are computationally 
expensive to evaluate. We are thus bound to small or intermediate chemical systems. 
The expansion is nonexistent in solids. 

(b) The single-determinant Slater-Jastrow wave function $> T ~ J provides in many cases 
good results for energy differences such as cohesion and binding due to large can- 
cellation of errors. However, any single-determinant ^^T J , with exception of fully- 
polarized state, divides space into 4 (2 for each spin) nodal cells, while in many cases 
the exact ground state wave function has only two nodal cells (for more on conjecture 
of minimal nodal cell division, see Ch. [3]). 

(c) The accuracy of Slater determinant based wave function is strongly dependent on the 
quality of one-particle orbitals. These come from methods which do not incorporate 
the explicit correlation such as Jastrow correlation factor and therefore are usually not 
optimal. However, to optimize many orbital coefficients is still a formidable challenge. 

The above limitations forced us and many others to think about novel wave functions 
with capabilities beyond fy^T J . Since determinant is the simplest antisymmetric object 
constructed from one particle orbitals, what would this object look like if we would build 
it from two, three or more particle orbitals? For two particle orbitals this object is called 
Pfaffian and the resulting wave function a Pfaffian pairing wave function. 



correlation factor. The resulting Slater-Jastrow wave function ty. 



S-J 

T 




exp[£/ CO yf] 



Similarly to multi-determinantal case, we can construct the Pfaffian pairing wave 



function as a linear combination of several Pfaffians 




(2.39) 



Each ith Pfaffian in the expansion (2.39) 



Pfi{l, . . . , N) = A[4>i(h 2), 04(3,4), . . . , MN - 1, AO] (2-40) 
= pf [&(1,2), &(3,4), . . . , UN - 1, AO] (2-41) 
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is an antisymmetrized product of one type of N/2 pairing orbitals 0j. The description 
of Pfafnan functional form together with application of Pfaffian pairing wave functions to 
several systems and possible extensions is discussed in great detail in Ch. [4] 

Lets us go back to a generalization mentioned at the beginning of this section. For 
several decades, it was known that it is possible to improve the nodal accuracy of trial wave 
functions for homogeneous systems by introducing a quasi-particle coordinates X [TBI El 

Eg ng EH ED 123, i.e., 

R^X, (2.42) 

where Tbf represents the backflow transformation. As a consequence, the nodes of ^^(X) 
will be in general different from those of ^^(R). Recently, some progress was also reported 
with chemical (inhomogeneous) systems |23|, |24"] . In this dissertation we put these new 
developments into test on several systems. The quasi-coordinate of zth electron at position 
Ti is given as 

Xj = Ti +&(R) 

= r, + £f»(R) + £f(R) + Cf n (R), (2.43) 

where we have again divided the contributions to one-body (electron- nucleus), two-body 
(electron-electron) and three-body (electron-electron-nucleus) backflow terms. They can be 
further expressed as 

er(R.) = £xM r f (2.44) 
i 

tr(R) = Y /U (r ij )r ij (2.45) 

£f en (R) = ^^[wiirij^i^rj^Tij + w 2 (ri j ,r iI ,r jI )r iI ], (2.46) 

where r,,- = r j — r,- and %, u and w\ with W2 are similar to one, two and three-body Jastrow 
terms. The calculational and implementation details together with some results are further 
discussed in Ch. HJ 

2.5 Pseudopotentials 

Many properties of interest in electronic structure methods are well determined by 
the behavior of only the valence electrons. On the other hand, the computational cost for 
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QMC methods increases very rapidly with the atomic number Z (oc Z 5 ' 5 |70j to Z 6,5 [71J). 
It is therefore physically desirable and computationally necessary to approximate the effect 
of core electrons by pseudopotentials. 

The presence of core electrons of atoms complicates the calculation for two reasons. 
The shorter length scales associated with variation of wave functions near the atomic core 
with large Z require the decrease of a time step in QMC simulations. The other problem is 
the increase of fluctuations of local energy as Z gets bigger caused by improper cancellation 
of divergent kinetic and potential terms. 

Analogously to independent particle theories such as Hartree-Fock and DFT, in 
QMC we remove the core electrons of atoms from calculation by introducing effective core 
potentials, also called pseudopotentials. The action of pseudopotential is typically different 
for electrons with different angular momenta. It is conventional to divide the pseudopoten- 
tial operator to local (common to all angular momenta) and nonlocal (independent for each 
angular momentum) parts as 

V ps (K) = Vl oc (K) + V ntoc 

= + E*&M (2.47) 

i i 

where 

iff«(*0 = E*f(*OIO<i| (2-48) 

I 

sums the contributions to nonlocal pseudopotential over angular momenta I for ith electron. 
The portion of the local energy El(R) from action of V n i oc on ^t(R) is then 

# T (R) ^ # r (R) ' 

Each ith argument in the above sum is evaluated as angular integration over the surface of 
the sphere passing through the ith electron and centered on the origin. If we choose the 
z-axis parallel to r^, the contribution for the ith electrons gets 

^« - ETh)5±i x §te|^> ,, (,50 

where Pi denotes the Ith Legendre polynomial. The integral is evaluated numerically on a 
spherical grid of Gaussian quadrature providing exact results up to certain maximum value 
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of Imax ■ When the orientation of axes of a given quadrature are chosen randomly, the result 



is an unbiased Monte Carlo estimation of integral (2.50). The quadratures implemented in 
QWALK [3] use l max = 3 (Octahedron - 6 points on the sphere) or higher precision l max = 5 
(Icosahedron - 12 points). For more details, see original Ref. |72j . 

In the DMC method, the application of nonlocal operator V n \ oc causes the ap- 
proximate propagator (R| exp[— rV^; oc ]|R') to have a fluctuating sign. Consequently, the 
imaginary time Schrddinger equation with importance function /(R, r) = ^ / r(R)*I'(R, r) 
will have modified form 

<l/(R - 7) ]v 2 f(R,r) + V • [v D (R)/(R,r)] + [E L (R) - E T )f(R,r) 



dr 2 

fvwMR) VW oc #(R,r) 



(2.51) 



1 * T (R) *(R,r) J' 

where the last term is unknown. We can therefore introduce the localization approximation 



by neglecting the last term in Eq. (2.51). It has been shown by Mitas et. al, |72j, that 
the error in energy from localization approximation oc (^t — &o) 2 - Fortunately, accurate 
Jastrow-based trial wave functions are available in most cases, so the localization error is 
small. Comparison of calculations with experiments for small systems had shown that the 
error from localization approximation is typically smaller then from fixed-node approxima- 
tion. Further tests on transition metal atoms |73[ 174] revealed some dependence of errors 
from localization approximation on quality of trial wave functions. Recently, a new algo- 
rithm based on lattice-regularization by Causula [75] was demonstrated to eliminate the 
localization approximation. 



2.6 Optimization of Variational Wave Functions 

As we have argued earlier in this chapter, the quality of trial wave functions 
controls the efficiency of VMC and DMC methods as well as final accuracy of their results. 
Given the freedom of choice for a wave function form in Monte Carlo, the accurate wave 
function may depend on many linear and non-linear parameters. The challenge is then the 
effective optimization of these variational parameters. 
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2.6.1 Variance Minimization 

Let us denote the set of variational parameters {c} of some real valued trial vari- 
ational wave function ^{ c }- The mean value of local energy (or VMC energy) with respect 
to ^{ c } is given as 

4?} = J < c > = IE c >) = (2.52) 

J {c} 

The variance of local energy is given as 



a 



2 J^ c} (4 c> -^ {c } ) 2 dR _ /fjp{c} 



((E^-E^) 2 ), (2.53) 



where we introduced a notation of the form: 



fO(R)^ 2 dR 

The simplest of the wave function optimization methods is to minimize the vari- 



ance [69 j in Eq. (2.53) on the set of fixed, finite MC configurations, where the walkers are 



distributed according to ^ 2 C0 y The set of starting variational parameters is denoted as 
{co}. The variance is generally a good function to minimize, since it is always positive and 
bounded from bellow by zero. In practice, the mean value of local energy at the end of 



optimization is not a priori known, so we replace in Eq. (2.53) by a reference energy 



E ref . 

There are several modifications to variance minimization. One standard improve- 
ment is to use weights, which account for wave function change after each step of opti- 
mization. The re-w heightened variance with the new set of parameters {c new } is then given 
as 

/ ^{c }^{^^i( R )(4 C " em} - Eref? dR 

where we have introduced weights 

*? x( R ) 

W (c )(R) = :T ' . (2-56) 

The advantage of the re-weighting scheme is a more accurate value of variance at each step 
of optimization. However, for optimization on small MC samples, sooner or later one weight 
will start to dominate over the others thus biasing the variance estimate. 
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Given the tools above, it is now possible to employ a suitable minimization algo- 
rithm for the search of a minimum. We use two different methods. First one is a modified 
version of a quasi- Newton method [76 . It uses only the value of a function to be minimized 
and builds up the information about the curvature of parameter space throughout the min- 
imization. Second algorithm is based on the general Levenberg-Marquardt (LM) method 



discussed bellow (see Sec. 2.6.3). It has a smooth transition between Newton method and 



steepest descent method with build-in stabilization parameter. However, it requires the in- 
formation about gradient and Hessian of a minimized function and therefore it is more costly 
to calculate. The gradient of variance with respect to ith. parameter is readily obtained by 



differentiating Eq. (|2.53|) as 

a a = 2 



(E L ,{E L -E)) + {-^E 



El 



2E 



(E L - E) 



(2.57) 



where subscript i denotes Since the variance minimization method can be viewed as 
a fit of the local energy on a fixed MC configurations [HS], an alternative expression for 
variance gradient follows from ignoring the change of the wave function 



a 2 i =2{E L AE L -E)). 



(2.58) 



The Hessian derived from gradient (2.58) is then given as 



af j = 2((E L)i -E)(E LJ -E)). 



(2.59) 



Important characteristic of Hessian (2.59) is that it is symmetric in i and j and positive 
definite. 



2.6.2 Energy Minimization 



A straightforward minimization of mean energy Eq. (2.52) is in general quite un- 
stable. The reason is that for sufficiently flexible variational wave function it is possible to 
lower the minimum of energy for a finite set of MC configurations, while in fact raising the 
true expectation value. As we have mentioned earlier, the variance is bounded from bellow 
and therefore this problem is far less severe |69j . The simple implementation therefore re- 
quires to minimize the energy on large MC configurations and possibly adjust for parameter 
changes by re-weighting procedure. Even then, as the initial and final wave function start 
to diverge, it is necessary to re-sample the MC configurations. 
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Similarly to variance minimization, we can improve the convergence of the op- 
timization by employing LM method and use information about gradient and Hessian of 
mean of local energy. The gradient of the mean of local energy can be readily obtained 



from Eq. (2.52 ) as 



E; 



—E L + 



2E- 



E)) 



(by Hermicity). 



(2.60) 



Note, that the expression in the last step of Eq. (2.60) has a favorable property of zero 
fluctuations as — > <l>o- Taking direct derivative of Eq. (2.60), the Hessian is 



E, 



+ 



E: 



(E L 



E) 



El, 



(2.61) 



* / * \ * 

It is clear that the above Hessian is not symmetric in i and j when approximated by 
finite sample. However, Umrigar and Filippi [77] recently demonstrated that the fully 
symmetric Hessian written entirely in the terms of covariances ((ab) — (a)(6)) has much 



smaller fluctuations then Hessian (2.61). Following their approach, modified Hessian from 
Ref. [77] is then given as 



E, 



+ 



* 2 



\ * 



{E L - E) 
(El,]) 



+ 



(Hi 

\ * 



El 



\ * 



E; 



(E L ,i) 



(2.62) 



L\ * ' 7 \ * 

where we added three additional terms (^-EL t i\ (^Elj) and (-|f£x,j) of zero expectation 
value (for proof, see e.g. Ref. [78 J. The first term makes the Hessian symmetric, while the 
remaining two terms put Hessian into covariance form. Hence, the addition of terms of zero 
expectation value for infinite sample has the effect of cancellation of most of the fluctuations 
in finite sample making the minimization method vastly more efficient. 



Another useful rearrangement of the Hessian (2.62) is 



E, 



'j 



+ 2 



+ 



r/*< 



* 2 



(E L - E) 



^1 



El, 



+ 



(El 
\ ^ 



\ ^ 



(El, 



(2.63) 
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If our minimization procedure involves only the linear parameters of the form exp[— c(i)f] 



(e.g. linear coefficients in Jastrow factor), the two terms on the first line of Eq. (2.63) cancel 
out, thus removing the need for expensive calculation of the Hessian of a wave function. 

In analogy to variance minimization, once the gradient and Hessian are computed, 
we search for a new set of parameters with lower value of energy using the Levenberg- 



Marquardt method discussed bellow (see Sec. 2.6.3). 



2.6.3 Levenberg-Marquardt Method 

Levenberg [79] and Marquardt [80] suggested to use a damped Newton method 
where the new step Ac in the parameter space is given by 

(H + M I)Ac = -g, (2.64) 

where H and I are Hessian and identity matrices, g denotes a gradient and \i is the positive 
damping parameter. This scheme has several favorable properties. 

(a) Positive \x makes H+/iI positive definite, which ensures that Ac is in descent direction. 

(b) For large values of [i we get 

Ac~--, (2.65) 
A* 

which is a short step in steepest decent direction. This approach is favorable, if we 
are far from minimum and the Newton approximation is not valid. 

(c) If fi is very small, then the new step Ac is in the Newton direction, which if we are 
close to minimum leads to (almost) quadratic convergence. 

What is left to determine is the optimal value for the damping parameter fx. The 
choice of \i is related to the values of Hessian H. If the Hessian is not positive definite (which 
can happen for small samples or if we are far from minimum) , we add to Hessian diagonal 
an absolute value of its lowest negative eigenvalue. Further, we start with /io = rmax(Hjj), 
where the choice of r is determined by the distance from minimum (from r = 1CP 6 if we 
are very close, up to even r = 1). Then we perform a correlated MC run (VMC or DMC) 
with three different damping parameters (e.g. ji\ = fiQ, [i<i = lO/^o and ^3 = 100/io)- A 
new damping fi new is then chosen at the point where the parabola fitted to the three new 
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energy values attains minimum (p new £ {pi, ££3)). On the other hand, if the correlated MC 
run is too expensive, we adjust the p based on the value of the gain ratio 

±Ac(Ac - g) 

where F(c) is the function of parameters c to be minimized. The value of a new damping 
p new is determined in the following way. 

• If p > 

1 , 
Pnew = iiomaxf-, 1 - (2p - 1) ], (2.67) 

• else 

Pnew = 2 n //0) (2.68) 

where n denotes a number of consequent times p < 0. For more details on Levenberg- 
Marquardt method and other least squares methods see for example Ref. [81. 



2.6.4 Implementation of Optimization Methods in QWALK 

Historically, the first implementation of minimization method in QWALK code 
was the OPTIMIZE method, which uses only the value of a function and a quasi-Newton 
minimizer. The optimization is bound to a fixed MC sample. It has proved effective for 
variance minimization of linear and nonlinear Jastrow parameters, however, the energy 
minimization required large MC samples. 

This deficiency was partially removed, when we have added the OPTIMIZE2 
method, which uses modified Hessians and gradients of variance and energy of Ref. |77j 
together with the Levenberg-Marquardt minimization algorithm. While the choice of mod- 
ified Hessian and gradient decreases significantly their fluctuations, the additional informa- 
tion (they provide) consequently improves the convergence of the method. The down-side 
to this approach is an additional cost of their calculation. For this reason, the analytical 
gradients and Hessians of wave function with respect to selected variational parameters 
were implemented (determinantal weights and orbital coefficients, Pfaffian weights and pair 
orbital coefficients). OPTIMIZE2 performs on a fixed MC sample, so the damping pa- 
rameter is determined via the gain ratio p as described in Levenberg-Marquardt method, 
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Table 2.1: List of some useful cost functions for a wave function optimization with denoted 
implementation in QWALK [3 . 



Cost Function 


Minimized quantity QWALK implementation 

OPTIMIZE OPTIMIZE2 NEWTON_OPT 


Variance 


m - m 


V 


V 


V 


Energy 


E 


V 


V 


V 


Mixed 


xE+{l-x)((E L -E) 2 ) 


V 


V 


V 


Absolute value 


(\El~E\) 


V 






Lorentz 


(ln(l + (E L -E) 2 /2)} 


V 






Ratio 


((E L -m 

E 








Overlap 


J * t *(t^oo) 









Sec. 2.6.3 OPTIMIZE2 has proved to be very effective method and was used for most of 
the optimizations in this dissertation and in other references [82 . 

As an attempt to remove the dependence of OPTIMIZE2 method on the fixed 
MC sample, I have also implemented the NEWTON.OPT method. Its principal advantage 
is that the energy and variance are obtained from a small MC run performed after each 
step of optimization. In addition, the optimal dumping can be chosen by an additional 



correlated MC run as described in Levenberg-Marquardt method, Sec. 2.6.3 Increased 
stability, however, comes at the price of additional MC calculations. 

Besides the well-known variance and energy minimization, it is possible to minimize 
several other so-called cost functions. Their incomplete list can be found in the Table |2.1| 
As was pointed recently by Umrigar and Filippi |77j . the wave functions optimized by 
minimization of a mixture of energy and variance have almost as low energy as energy 
optimized wave functions but a lower variance. 

We demonstrate this for the ground state wave function of N2 molecule (see 



Fig 2.2). The employed NEWTON.OPT method minimized 23 Jastrow parameters (the 
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curvatures of basis functions and all linear coefficients) of single determinant Slater-Jastrow 
wave function on 56000 walkers. The mixture minimized energy is almost as good as energy 
minimized energy, while the mixture minimized dispersion a falls in between a from energy 
minimization and variance minimization. Also note rather slower converge of energy and 
mixture minimization when compared to variance minimization. This behavior may have 
two reasons. One is the more complex structure of energy function. Alternatively, we just 
happen to be much further from minimum for energy function than for variance. 

The mixture minimized wave functions are therefore the most efficient for DMC 
calculations. In the same time, small addition of variance (typically 5%) to energy in 
minimization greatly decreases the fluctuations of Hessian and enables us to use smaller 
MC samples. 

2.7 Summary 

In this chapter, we tried to plot an overview of two principal Quantum Monte Carlo 
methods, the variational and diffusion Monte Carlo. Further, we mentioned some essential 
sampling MC techniques such as Metropolis, importance and correlated sampling. A more 
technical details of QMC, like the use of non-local pseudopotentials are discussed. The only 
two uncontrolled approximations of DMC, the fixed node and localization approximations, 
are introduced to circumvent the fermion sign problem. We analyzed some of the properties 
of trial wave functions and presented two very general and accurate wave functions, the 
Slater-Jastrow and the Pfaffian-Jastrow trial wave functions. Both can be teamed-up with 
backflow transformation which can further lead to improved nodes. Finally, we talked about 
the optimization of variational trial wave functions and acquainted the reader with its three 
implementations in QWALK code. 
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Iteration number 




10 20 30 40 



Iteration number 

Figure 2.2: Energy E and dispersion of local energy a of N2 molecule versus iteration 
number of minimization. Notation: minimization of variance (green triangles), energy 
(blue circles) and 95% mixture of energy and 5% of variance (red squares). Upper figure: 
Energy versus iteration number. The error-bars are of the size of symbol or smaller. Inset: 
the later iterations on expanded scale. The mixture minimized energy is almost as good as 
energy minimized one while variance minimized energy is higher by almost 10 mH. Lower 
figure: Same as the upper figure but for the dispersion of local energy a rather than energy. 
Inset shows that the mixture minimized a falls in between a from energy minimization and 
variance minimization. 
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Chapter 3 



Nodal Properties of Fermionic 



Wave Functions 



Sections of this chapter also appeared in: 

Approximate and exact nodes of fermionic wave functions: 
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M. Bajdich, L. Mitas, G. Drobny, and L. K. Wagner, 
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In Condensed Matter Theories, vol. 20 (Nova Science Publishers, 2006). 
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3.1 Introduction 

Nodes of fermionic wave functions and their related objects, density matrices, are 
of great interest to physicist for several reasons. From QMC point of view, the knowledge 
of exact nodes of fermionic wave functions would enable us to obtain exact ground state 
energies by means of fixed-node DMC. Analogously, the knowledge of exact nodes of tem- 
perature density matrices would lead to the solution of the thermal fermion problem^in path 
integral Monte Carlo (PIMC). The nodal properties of wave functions dramatically change 
with increased correlations among electrons and are relevant to transport and many-body 
phases such as superconductivity. 

The general properties of fermion nodes were first analyzed in an extensive study 
by Ceperley [83_ , which included a proof of the tiling property and generalizations of the 
fermion nodes to density matrices. In addition, for some free particle systems, it was 
numerically shown [83] that there are only two nodal cells. The fermion nodes for degenerate 
and excited states were further studied by Foulkes and co-workers [S3]. Recently, Mitas 
used the property of connectivity from Ref. [83J to show that a number of spin-polarized 
noninteracting and mean-field systems (homogeneous electron gas, atomic states, fermions 
in the box, in the harmonic well and on the sphere) has ground state wave functions (given as 
Slater determinants) with minimal number of nodal cells [85, 86J. Further, he demonstrated 
that for spin-unpolarized systems an arbitrarily weak interaction introduced by Bardeen- 
Cooper-Schrieffer (BCS) wave function reduces the four non-interacting nodal cells to just 
two. Finally, he has also shown that the minimal number of nodal cells property extends 
to the temperature density matrices. 

The fermion nodes of small systems, mostly atoms, were investigated in several 
previously published papers [29] [55] [56] ET] [88] [89] [90]. Interesting work by Bressanini, 
Reynolds and Ceperley revealed differences in the nodal surface topology between Hartree- 
Fock and correlated wave functions for the Be atom explaining the large impact of the 
2s, 2p near-degeneracy on the fixed-node DMC energy [91]. More recently, improvement in 
fixed-node DMC energies of small systems using CI expansions |92[ [93] , and pairing wave 
functions [HI \TE\ [M] [§S] were also reported. 



This chapter is organized as follows. In Sec. 3.2 we summarize the general prop- 



1 The thermal fermion problem is the construction of a polynomial-time algorithm of the exact thermody- 
namics of many- fermion system at positive temperature |83| . 
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erties of fermion nodes. In Sec. 3.3 we discover new exact fermion nodes for two and 



three-electron spin polarized systems. In Sec. 3.4 we categorize the nodal surfaces for the 



several half-filled subshells relevant for atomic and molecular states. In Sec. 3.5 we show 
how opposite spin correlations eliminate a redundant nodal structure of HF wave functions 
for two specific cases of spin-unpolarized states. Finally, in the last section we present our 
conclusions and suggestions for future work. 



3.2 General Properties 

Let us assume a system of spin-polarized electrons described by a real wave function 
VP(R). Then the ^(R) is antisymmetric with respect to any particle exchange as described 
in Sec. |2.4[ Consequently, there exists a subset of electron configurations {R noc ; e }, called a 
fermion nodej^J for which the wave function is zero, i.e., 

^(Rnode) = 0. (3.1) 

Naturally, we eliminate from the definition the regions, where the wave function vanishes be- 
cause of other reasons (e.g., external potential) than antisymmetry. In general, the fermion 



node is a (ND — l)-dimensional manifold (hypersurface) defined by implicit Eq. (3.1) as- 
suming that we have N fermions in a D— dimensional space. When the positions of any two 
electrons are equal (i.e., r,j = rj), the antisymmetry ensures that the ^(R) = 0. However, 
that does not fully specify the nodes, but only defines the (ND — D)-dimensional subspace 
of coincidence planes, where wave function vanishes. Note that when we talk about the 
nodes, we always mean the nodes of a many-body wave functions and not the nodes of 
one-particle orbitals. 

Let us now introduce the basic properties of fermion nodes as they were studied 
by Ceperlev[83] some time ago. 

(a) Tiling property for the non- degenerate ground state — Let us define a nodal cell O(R^) 
as a subset of configurations, which can be reached from the point R^ by a continuous 
path without crossing the node. The tiling property says that by applying all possible 
particle permutations to an arbitrary nodal cell of a ground state wave function one 



2 The first question about the node structure of fermion wave functions was raised by J. B. Anderson in 
1976 [56]. 
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covers the complete configuration space 91 (i.e., Yp fi(-PRt) + fi(Rt) = 9t). Note that 
this does not specify how many nodal cells are there. 

(b) Nodal crossings — If two nodal surfaces cross each other, they are orthogonal at the 
crossing. If n nodal surfaces cross each other, the crossing angles are all equal to ir/n. 

(c) Symmetry of the node — Symmetry of the state is also symmetry of the node. 

(d) Connectivity <^=^ two maximal nodal cells — It is possible to show that there are 
only two nodal cells using the argument based on connectivity of particles by triple 
exchanges. The three particles i, j, k are called connected if there exist a cyclic 
exchange path i — > j, j — ► k, k — > i, which does not cross the node. In addition, if all 
particles of some point Rt are connected together by triple exchanges, the VP(R) has 
only two nodal cells. In other words, the whole configuration space is covered by only 
one positive and one negative nodal cell (i.e., the nodal cells are maximal). This can 
be better understood if we realize the following. First, any triple exchange is just a 
two pair positive permutation. The existence of point R^ with all particles connected 
is however equivalent to have all positive (not flipping the sign) permutations P + R< 
belonging to the same nodal cell. Second, the tilling property implies, that once all 
particles are connected for Rt it is true for entire cell O(R^). Therefore, there will be 
only one maximal cell per each sign. More details on this property can be found in 
Ref. |83]. 

3.3 Exact Nodal Surfaces 

We assume the usual electron-ion Hamiltonian and we first investigate a few- 
electron ions focusing on fermion nodes for subshells of one-particle states with s,p,d, /... 
symmetries using variable transformations, symmetry operations and explicit expressions 
for the nodes. 

3.3.1 Three-Electron Quartet 4 S(p 3 ) State 

Let us first analyze a special case with n = ri and r<iz = r^\. It is then easy 
to see that the inversion around the origin with subsequent rotations is equivalent to the 
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Figure 3.1: Inversion and two subsequent rotations of three particles: (a) Original and 
inverted (primed) positions; (b) Positions after the rotation R\ in the plane given by the 
particles 1, 2 and the origin; (c) Positions after the second rotation R2 around the ri + r 2 
axis. Note that the original positions of the particles 1 and 2 are exchanged. 



exchange of two particles, say, 1 and 2 (Fig. 3.1 ). Therefore, for this particular configuration 
of particles, the combination of parity and rotations is closely related to the exchange 
symmetry. The illustration also shows that the six distances do not specify the relative 
positions of the three electrons unambiguously. For a given set of the distances there are 
two distinct positions, say, of the electron 3, relative to the fixed positions of electrons 1 



and 2 (see Fig. 3.1 ) and compare positions 3 and 3" of the third electron. 



In order to analyze the wave function in an unambiguous manner it is convenient 



to define new coordinates. Let us denote r 



12 



ri + r 2 ,r 



12 



r^ 2 |, together with the 



customary ri2 = ri — r2,ri2 
Cartesian coordinates 



I r 12 1 • We can now introduce the following map of the 



(ri,r 2 ,r 3 ) -» (rf 2 , n 2 , r 3 , cos a, cos /3, 7, f2) 



(3.2) 
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with definitions: cosa = r3 • (n x r2)/(r3|ri x r2 1) , cos/? = rf 2 • ri2/(^i2 r i2) arid 7 being 
an azimuthal angle of r3 in the relative coordinate system with unit vectors = r^/r^, 
e z = r i x r2/|ri x r2 1 , e y = e z x e x . For completeness, O denotes three Euler angles, which 
fix the orientation of the three-particle system in the original coordinates (e.g., two spherical 
angles of ri x r2 and an azimuthal angle of r^ 2 ). Since the angles f2 are irrelevant in S 
symmetry, the first six variables fully specify the relative positions of the three particles and 
the wave function dependence simplifies to ^(rf 2 , t*12 ; r 3, cos a, cos /?, 7). Consider now two 
symmetry operations which change the sign of the wave function and keep the distances 
unchanged: parity Pj and exchange P12 between particles 1 and 2. The exchange flips the 
sign of all three cosa, cos (3, 7 while the parity changes only the sign of cosa. The action 
of P/P12 on ^ leads to 

cos a, — cos /?, —7) = cos a, cos /?, 7) (3-3) 

showing that the wave function is even in the simultaneous sign flip (cos /3, 7) — > (— cos (3, —7). 
Applying the exchange operator P12 to the wave function and taking advantage of the pre- 
vious property gives us 

\&(..., — cos a, cos /?, 7) = — ^(..., cos a, cos (3, 7) (3-4) 

suggesting that there is a node determined by the condition cos a = 0. It is also clear that 
the same arguments can be repeated with exchanged particle labels 2 *-* 3 and 3 *-* 1 and 
we end up with the same nodal condition, r% ■ (17 x T2) = 0. This shows that the node is 
encountered when all three electrons lie on a plane passing through the origin. Now we 
need to prove that this is the only node since there might possibly be other nodal surfaces 
not revealed by the parametrization above. The node given above clearly fulfills the tiling 
property and all symmetries of the state. Furthermore, the state is the lowest quartet of 
S symmetry and odd parity (lower quartets such as Is2s3s,ls2s2p, and ls2p 2 have either 
different parity or symmetry) and for the ground state we expect that the number of nodal 
cells will be minimal. This is indeed true since the node above specifies only two nodal cells 
(one positive, one negative): an electron is either on one or the other side of the nodal plane 
passing through the remaining two electrons. Furthermore, any distortion of the node from 



the plane necessarily leads to additional nodal cells (see Fig. 3.2), which can only increase 
energy by imposing higher curvature (kinetic energy) on the wave function. This is basically 
the Feynman's argument from the proof demonstrating that the energy of fermionic ground 
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Figure 3.2: (a) An illustration of an artificial distortion of the planar ground state node 
for the quartet state (see text); (b) The original and parity transformed distorted node; (c) 
Finally, a subsequent rotation of the inverted distortion necessarily leads to a new nodal 
pocket which is artificial for the ground state. In fact, nodes with similar topologies are 
present in excited states (see bellow). 

state is always above the energy of the bosonic ground state (and also essentially the same 
argument as used for the proof of the tiling property [83J. In fact, all higher excited states 
of this symmetry (HF wave function of 2p orbitals, e.g., excited state of B +2 ion) have 



additional nodes, as expected (see Fig. 3.3). Given all the arguments above we conclude the 
proof that the plane is the exact node. Note that it is identical to the node of Hartree-Fock 
wave function of 2p orbitals given by ^>hf = det[p(r)x, p(r)y, p(r)z] where p[f) is a radial 
function. 

The coordinate transformation above is not the only one that can be used to 
analyze this state. The high symmetry of the problem enables us to find an alternative 
coordinate map with definitions of cos modified to cos /?' = [(ri x Y2) x r^ 2 ] -r^/fKri xr-2) x 
r 12 1 1 r i2 1] an d 7 to 7' by redefinition of e z to e' z = [[(t\ x Y2) xri2] xr^~ 2 ]/| [(ri xr2) x 1*12] xr+ 2 | 
and e' y = e' z x e x . In the redefined coordinates the search for the node simplifies to an 
action of P\% on *&(rf 2 , f*i2> f3, cos a, cos 13', 7') 

$(..., — cos a, cos f3 , 7 ) = — ^(..., cos a, cos /?', 7') (3-5) 



since the distances and cos/?', 7' are invariant to P\2- Obviously, this leads to the same 
nodal condition as derived above. 
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Figure 3.3: The 3D projected nodes of a few selected excitations for the symmetry adopted 
CI expansion of the 4 S'(p 3 ) ground state. The exact planar node of the quartet ground 
state is also possessed by all the excitations. The small spheres show the fixed positions of 
two electrons while the third one is scanning the nodal surface. Labels indicate the orbitals 
involved. 

It is quite interesting to compare these two coordinate maps with /?, 7 and /3',7'. 
Although parity and exchange are independent operators, the analysis above shows that in 
an appropriate coordinate system they imply the same nodal surface. Both these operators 
cause an identical sign change of the wave function indicating thus a special symmetry of the 
4 S'(p 3 ) ground state node, which is higher than would be expected solely from antisymmetry. 
Similar observation was made in a study of fermion node in another case of two electron 
atomic state [89j ES] . 
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3.3.2 Two-Electron Triplet 3 P(p 2 ) and 3 S 9 (vr 2 ) States 

Apparently, the exact node of this case was derived in a different context by Breit in 
1930 [88, [92], ES] • Here we offer an independent proof which enables us to apply the analysis 
to some molecular states with the same symmetries. The exact node for the 3 P(p 2 ) state 
can be found in a similar way as in the case of quartet above. The state has even parity, 
cylindric symmetry, say, around z-axis, and is odd under rotation by tt around x, y axes, 
R(ttx), R(iry). The mapping of Cartesian coordinates which enables to analyze the wave 
function symmetries is given by 

(ri,r 2 ) -> (rf 2 , ri2, cos u, cos [3, <p,(p') (3.6) 

where cos a; = zq ■ (ri x r2)/|ri x r 2 | with zq being the unit vector in the z-direction and 
ip' being the azimuthal angle of ri x r 2 ; <p' can be omitted due to the cylindric symmetry. 
Further, ip is the azimuthal angle of rf 2 in the relative coordinate system with the x-axis 
unit vector given by a projection of zo into the plane defined by ri,r2, i.e., e x = zo p /|zo p |, 
e z = (ri x r 2 )/|ri x r 2 | and e y = e z x e x . Action of P/Pi 2 i?(7rx) reveals that the wave 
function is invariant in the simultaneous change (cos/?, ip) — > (— cos 0, — p). This property 
and action of Pi 2 to the wave function together lead to 

#(...,- cos ...) = -^(...,coslj, ...) (3.7) 

with the rest of the variables unchanged. The node is therefore given by cos lo = and 
is encountered when an electron hits the plane which contains the z-axis and the other 
electron. As in the previous case the nodal plane fulfills the tilling property and manifestly 
divides the space into two nodal cells so that we can conclude that this node is exact. The 
exact node again agrees with the node of Hartree-Fock wave function = det[p(r)x, p{r)y\. 

The fixed-node QMC energies for the 4 S , (p 3 ) and s P(p 2 ) cases derived above were 
calculated for a nitrogen cation with valence electrons in these states. The core electrons 
were eliminated by pseudopotential [I] . The trial wave function was of the commonly used 
form with single HF determinant times a Jastrow correlation factor [9]. Note that the 
pseudopotential nonlocal s— channel does not couple to either odd parity S state or even 
parity P(p 2 ) state so that that the nonlocal contribution to the energy vanishes exactly. 

In order to compare the fixed-node QMC calculations with an independent method 
we have carried out also configuration interaction calculations with ccpV6Z basis [97] (with 
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Table 3.1: Total energies (in Hartrees) of N + , N +2 and N +3 ions with core electrons elimi- 
nated by pseudopotentials [lj. The energies are calculated by variational (VMC) and fixed- 
node diffusion (DMC) quantum Monte Carlo and configuration interaction (CI) methods. 
The HF energies are given as a reference for estimation of the correlation energies. 

State HF CI VMC DMC 

3 P(p 2 ) -5.58528 -5.59491 -5.59491(2) -5.59496(3) 

4 S"(p 3 ) -7.24716 -7.27566 -7.27577(1) -7.27583(2) 

5 S(sp 3 ) -8.98570 -9.02027 -9.01819(4) -9.01962(5) 



up to three g basis functions), which generates more than 100 virtual orbitals in total. In the 
CI method the wave function is expanded in excited determinants and we have included all 
single, double and triple excitations. Since the doubles and triples include two- and three- 
particle correlations exactly, the accuracy of the CI results is limited only by the size of the 
basis set. By comparison with other two- and three-electron CI calculations we estimate 
that the order of magnitude of the basis set CI bias is ~ 0.01 mH for two electrons and 
~ 0.1 mH. and for three electrons (despite the large number of virtuals the CI expansion 
converges relatively slowly [98] in the maximum angular momentum of the basis functions, 
in our case l max = 4). The pseudopotentials we used were identical in both QMC and CI 
calculations. 



The first two rows of Tab. 3.1 show the total energies of variational and fixed-node 
DMC calculations with the trial wave functions with HF nodes together with results from 
the CI calculations. For 3 P(p 2 ) the energies agree within a few hundredths of mH with 
the CI energy being slightly higher but within two standard deviations from the fixed-node 
QMC result. For 4 S(p 3 ) the CI energy is clearly above the fixed-node DMC by about 0.17 
mH as expected due to the limited basis set size. In order to illustrate the effect of the 
fixed-node approximation in the case when the HF node is not exact we have also included 
calculations for four electron state 5 S(sp 3 ) (for further discussion of this Hartree-Fock node 



see Sec 3.4.3 below). For this case, we estimate that the CI energy is above the exact one 
by ~ 0.3 mH so that the fixed-node energy is significantly higher than both CI and exact 
energies. Using these results we estimate that the fixed-node error is ~ 1 mH, i.e., close to 
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3% of the correlation energy. 

Since in the p 2 case we have assumed cylindric symmetry, the derived node equa- 
tion is applicable to any such potential, e.g., equidistant homo-nuclear dimer, trimer, etc, 
with one-particle orbitals vr x ,7r y which couple into the triplet state s T, g (Tr x 7r y ). 

Note that the parametrization given above automatically provides also one of the 
very few known exact nodes in atoms so far |29| [57] . i.e., the lowest triplet state of He 
3 S(ls2s). The spherical symmetry makes angles u> and tp irrelevant and simplifies the two- 
electron wave function dependence to distances r\,r2,r\2 or, alternatively, to r\2, rf 2 , cos/3. 
Applying P12 to wave function ^(ri 2 , rf 2 , cos (3) leads to 

- *(ri 2 ,rf 2 ,cos/3) = ^(r 12 ,rf 2 ,- cos /3) (3.8) 

so that the node is given by the condition cos [3 = 0, i.e., r\ — r 2 = 0. 

In addition, the presented analysis sheds some light on the He 3 P(ls2p) state 
node which was investigated before |89j as having higher symmetry than implied by the 
wave function symmetries. The symmetry operations reveal that the wave function de- 
pends on I cosu;| and that the node is related to the simultaneous flips such as (cos j3, (p) —* 
(— cos (3, —(f) or angle shifts ip —>(/? + tt. Since, however, two of the variables are involved, 
the node has a more complicated shape as the previous study illustrates [89 . In order to 
test the accuracy of the HF node we have carried out a fixed-node diffusion Monte Carlo 
calculation of the He 3 P(ls2p) state. The resulting total energy of -2.13320(4) H is in an 
excellent agreement with the estimated exact value of -2.13316 H [99 , which shows that 
the HF node is very close to the exact one [100 . 

3.4 Approximate Hartree-Fock Nodes 

It is quite instructive to investigate the nodes of half-filled subshells of one-particle 
states with higher angular momentum. 

3.4.1 Hartree-Fock Node of 6 S(d 5 ) State 

The HF determinant wave function for 6 S*((i 5 ) is given by 
5 

^hf = JJp(ri)det[2z 2 - x 2 - y 2 ,x 2 - y 2 ,xz,yz,xy], (3.9) 

i=l 
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Figure 3.4: The 3D projected Hartree-Fock node of 6 S(d 5 ) state, which is an elliptic 
cone (left and right pictures). The middle picture illustrates a case when two pairs of two 
electrons lie on orthogonal planes which pass through the origin. This two-plane node is of 
lower dimension because of the additional condition on positions of the electrons. It appears 
as a crossover between the cones with different orientation (left and right pictures). The 
small spheres show the positions of the four electrons while the line denotes the z— axis. 

where p{ri) is the radial part of the d-orbital and we assume that all the orbitals are from 
the same I = 2 subshell, e.g., 3d subshell. Since all radial functions are the same they factor 
out from the determinant and for the purpose of finding the node they can be omitted. 
The S symmetry allows to rotate the system so that, say, electron 1 is on the z-axis, and 
then the corresponding column in the Slater matrix becomes (2zf,0, 0, 0, 0). Assuming that 
z\ ^ we can then write the nodal condition as 

det[x 2 — y 2 , xz, yz, xy] = 0. (3.10) 

Using one of the electrons as a probe (i.e., looking at the node from the perspective of one 
of the electrons) we can find the projection of the node to 3D space. By denoting the probe 
electron coordinates simply as (x, y, z) and by expanding the determinant we get 

(x 2 — y 2 )m\ + xzm2 + yzms + xym^ = (3-11) 

where m; are the corresponding cofactors. We divide out the first cofactor assuming that 
it is nonzero (not a crucial assumption as clarified below). We get 

(x 2 — y 2 ) + axz + byz + cxy = (3-12) 

where o = TO2/7711, b = m^/mi, c = m^jm\. By completing the square this can be further 
rearranged to 

(x — k\y){x — kiy) + z{ax + by) = (3.13) 



3.4. Approximate Hartree-Fock Nodes 



45 



with ki j2 = (— c± Vc 2 + 4)/2. Let us define rotated and rescaled coordinates 

u* = -{ak 2 - b)(x - k iy )/{h - k 2 ) (3.14) 

v* = (ak 1 -b)(x-k 2 y)/(fa-k 2 ) (3.15) 

w* = z[{aki - b){ak 2 - b)}/{ki - k 2 f (3.16) 



so we can write the Eq. (3.12) as 

u*v* + w*u* + w*v* = 0. (3.17) 



Note that this equation has a form which is identical to Eq. (3.11) with mi = so this 



representation is correct for general mi. After some effort one finds that Eq. (3.17) is a 



cone equation (i.e., d z 2 orbital) as can be easily verified by using the following identity 

(2u 2 - v 2 - w 2 )/8 = u*v* + w*u* + w*v*, (3.18) 

where u = u* + v* + 2w* , v = (-u* + v* + 2w*), w = (u* - v* + 2w*). The 3D projected 
node is therefore rotated and rescaled (elliptic) cone. 

At this point it is useful to clarify how the derived node projection cone is related to 
the complete 14-dimensional node. Remarkably, the 3D projection enables us to understand 
some of the properties of the 14-dimensional manifold. First, the cone orientation and 
elliptic radii (i.e., rescaling of the two axes with respect to the third one) are determined by 
the position of the four electrons in 3D space: with the exception of special lower dimensional 



cases explained below there always exists a unique cone given by the Eq. (3.17) which "fits" 
the positions of the four electrons. Besides the special cases (below) we can therefore define 
a projection of a single point in 4 x 3 = 12-dimensional space of four electrons onto a cone. 
That also implies that the complete 12-dimensional space describes a set (or family) of cones 
which are 3D projections of the nodal manifold. Similar projection strategies are often used 
in algebraic geometry to classify or analyze surfaces with complicated topologies or with 
high dimensionality. 

Since the cone orientation and two radii are uniquely defined by the point in 12 
dimensions and the cone itself is a 2D surface in 3D space of the probe electron, the complete 
node then has 12+2=14 dimensions. Therefore the d 5 HF node is a set of cone surfaces 
specified by the positions of the electrons . This particular form is simply a property of 
the d 5 Hartree-Fock determinant. From the derivation above it is clear that after factoring 
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out the radial parts one obtains a homogeneous second-order polynomial in three variables 
with coefficients determined by the positions of the four electrons. In fact, from the theory 
of quadratic surfaces |1U1| . one finds that a general elliptic cone can possibly fit up to 
five 3D points/electrons, however, in our case the cone has an additional constraint. Our 
system was reoriented so that one of the electrons lies on the z-axis; that implies that the 
z-axis lies on the cone. Therefore the cone always cuts the xy (i.e., z = 0) plane in two 
lines, which are orthogonal to each other. The orthogonality can be verified by imposing 



in Eq. (3.13) and checking that k\k2 = — 1. In addition, one can find "degenerate" 



configurations with two pairs of two electrons lying on orthogonal planes (Fig. 3.4). This 
corresponds to the "opening" of the cone with one of the elliptic radii becoming infinite 



and the resulting node having a form of two orthogonal planes (Fig. 3.4). Since in this 
case there is an additional condition on the particle positions, the two-plane node has lower 
dimension and is a zero measure subnode of the general 14-dimensional node. The condition 
is equivalent to A44 = b 2 — a 2 — abc = 0, where A44 is one of the quadratic invariants [101 . 
There are more special cases of lower dimensional nodes: (a) when two electrons lie on a 
straight line going through the origin; (b) when three electrons lie on a plane going through 
the origin; (c) when four electrons lie in a single plane. 

Remarkably, the analysis above enables us to find the number of nodal cells. From 



Fig. 3.4 one can infer that by appropriate repositioning of the four electrons the cone 
surface smoothly "unwraps" the domains inside the cone, forms two crossing planes and 
then "wraps" around the cone domains of the opposite sign. That implies that an electron 
inside one of the cone regions can get to the region outside of the cone (with the same wave 
function sign) without any node crossing, using only appropriate concerted repositioning 
of the remaining four electrons. That enables us to understand that a point in the 15- 
dimensional space (positions of five electrons) can continuously scan the plus (or minus) 
domain of the wave function: there are only two maximal nodal cells. 



3.4.2 Hartree-Fock Nodes of the 8 S(f) Ion 

We will use similar strategy as in the preceding case. After rotating one of the 
electrons to z-axis we expand the determinant in the probe electron column and eliminate 
the radial orbitals which form an overall prefactor of the Slater determinant since we assume 
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cone times planar surface or a cone "fused" with planar surface what forms a single sheet 
surface. There is a smooth transition between these two forms depending on the positions 
of six electrons which are denoted by the small spheres. Note that the node contains the 
z— axis which is denoted by the dashed line. 

that all seven /-states are from the same I = 3 subshell (e.g., 4/). We get 

(mix + m2y)(4:Z 2 — x 2 — y 2 ) + m^zix 2 — y 2 ) + m^xyz 

+m 5 x(x 2 - 3y 2 ) + m 6 y(y 2 - 3x 2 ) = 0. (3.19) 

Note that the node contains the z-axis and there are two possible values of z for any x, y 
since the form is quadratic in z. This restricts the node shapes significantly and by further 
analysis one can find that the nodal surface projection into 3D has two topologies (Fig. |3.5j ). 
The first one is a cone times a planar surface (topologically equivalent to the I40 spherical 
harmonic). Note that, in general, the planar surface is deformed from a straight plane since 
it passes through the origin and, in addition, it fits three of the electrons. The second 
topology is a "fused" cone and planar surface, which results in a general single sheet cubic 
surface. The node transforms smoothly between these two topologies depending on how 
the six electrons move in space. These two topologies define the projection of the node into 
the probe 3D space and therefore enable us to capture the many-dimensional node for this 
particular Hartree-Fock state. This again enables to describe the complete node using a 
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Figure 3.6: The 3D projection of the nitrogen cation 5 S(sp' i ) Hartree-Fock node (the core 
electrons are eliminated by pseudopotentials). The projected node exhibits two topologies. 
It is either a planar surface deformed by the radial orbital functions at the nucleus or, 
in certain configurations, the deformation forms a small bubble detached form the surface 
(the picture on the right). The small cross is the location of the ion while the small spheres 
denote positions of electrons. 

theorem from algebraic geometry which states that any cubic surface is determined by an 
appropriate mapping of six points in a projective plane |102| I103[ 1104] . To use it we first need 
to realize the following property of the 3D projected node: The node equation above contains 
only a homogeneous polynomial in x,y,z which implies that in spherical coordinates the 
radius can be eliminated and the node is dependent only on angular variables. Hence, any 
line defined by an arbitrary point on the node and the origin (i.e., a ray) lies on the node. In 
other words, we see that the surface is ruled, i.e., it can be created by continuous sweep(s) 
of ray(s) passing through origin. This enables us to project the positions of the six electrons 
on an arbitrary plane, which does not contain the origin, and the node will cut such a plane 
in a cubic curve. As we mentioned above, a theorem from the algebraic geometry of cubic 
surfaces and curves says that any cubic surface is fully described by six points in a projective 
plane (see \1U'2\ I1L)3| I1U4] ) . For ruled surfaces any plane not passing through the origin is 
a projective plane and therefore we can specify a one to one correspondence between the 
6 x 3 = 18 dimensional space and our cubic surface in 3D. Obviously, there will be a number 
of lower-dimensional nodes which will correspond to positions of electrons with additional 
constraints such as when they lie on curve with the degree lower than cubic; i.e., a conic. 
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3.4.3 Hartree-Fock Nodes of the 5 S(sp 3 ) Ion 

The HF node for this two-shell spin-polarized state can be investigated in a similar 
way as in previous cases with a new feature that the radial parts will be present in the 
expansion of the determinant. By expanding the determinant in the column of the probe 
electron with position x, y, z the 3D node projection is simply given by 

x + b'y + c'z + d'r](r) = 0, (3.20) 

where b',c',d' depend on ratios of cofactors and r]{r) = p s {r) / p p {r) is the ratio of radial 
parts of s and p orbitals and r = \J x 2 + y 2 + z 2 . The probe electron will see a plane with a 



approximately bell-shape deformation in the area of the nucleus (see Fig. 3.6 ). The shape of 
deformation depends on the ratio of s and p radial parts and the magnitudes and signs of the 
cofactors. For certain configurations the deformation is so large that it gets detached from 
the surface and forms a separated ellipsoid-like bubble. The bubble results from the radial 
dependence of rj(r) which for pseudized core is not a monotonic function and therefore can 
create new topologies. Note that despite the fact that the 3D projection shows a separated 
region of space (the bubble) the complete node has again the minimal number of nodal cells 
property. To understand this, suppose that the probe electron is located inside the bubble 
and wave function there has a positive sign. Let us try to imagine how the electron can get 
to the other positive region (the other side of the planar surface). Seemingly, the electron 
would need to cross the nodal surface twice (the surface of the bubble and the planar 
surface). However, the complete node is a collective-coordinate object and by moving the 
other two electrons in an appropriate way the bubble attaches to the surface and then fuses 



into a single surface (Fig. 3.6, left) so that the probe electron can reach the positive region 
without node crossing. 

In order to see whether the correlation would change the HF node we have car- 
ried out a limited study of the CI wave function nodes for this case; we have found some 
differences but we have not discovered any dramatic changes to the HF nodes. To quantify 
this further we have calculated the CI energy (with the same basis and level of correlation 



as in the previous cases) and the result is in the last row of Tab. 3.1 We estimate that the 
fixed-node bias of the HF node is of the order of ~ 0.001 Hartree which is close to ~ 3% of 
the correlation energy. Obviously, the DMC energy is above the exact one and percentage- 
wise the amount of missing correlation energy is not insignificant. We conjecture that the 
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HF node is reasonably close to the exact one although the fine details of the nodal surface 
are not captured perfectly. 

3.4.4 Hartree-Fock Nodes of Spin-Polarized p 3 d 5 and sp 3 d 5 Shells with S 
Symmetry 

Let us for a moment assume a model wave function in which the radial parts of 
s,p, d orbitals are identical. Then, using the arrangements similar to d 5 case, we can expand 
the determinant of p 3 d 5 in one column and for the 3D node projection we then get 

2u 2 - v 2 - w 2 + au + 0v + = 0, (3.21) 

where u, v, w are appropriate linear combinations of x, y, z. This can be further rewritten 
as 

2(u + a/4) 2 - (v - P/2) 2 -(w- 7 /2) 2 + 5 = 0, (3.22) 

where 6o = (— a 2 /2 + f5 2 + 7 2 )/4. It is clear that the quadratic surface is offset from the 
origin (nucleus) by a vector normal to au + (3v + jw = plane. Using the properties of 
quadratic surfaces one finds that for (a 2 /(a 2 + /? 2 + 7 2 )) < 2/3 the node is a single-sheet 
hyperboloid with the radius Vooj otherwise it has a shape of a double-sheet hyperboloid. 
The double-sheet hyperboloid forms when there is an electron located close to the origin. 
A special case is a cone which corresponds to (<5q = 0). The case of sp 3 d 5 is similar, but 



with different 5q, which now has a contribution from the s-orbital (see Fig. 3.7). Once we 
include also the correct radial parts of orbitals in the s,p,d channels the coefficients of the 
quadratic form depend on both cofactors and orbital radial functions. The resulting nodal 
surface is deformed beyond an ideal quadric and shows some more complicated structure 



around the nucleus (see Fig. 3.8) as illustrated on HF nodes of the majority spin electrons 
in Mn +2 ion (note that the Ne-core electrons were eliminated by pseudopotentials). 



3.5 Nodes and Spin Correlations 

We conjecture that a non-degenerate ground state of any given symmetry possesses 
only two maximal nodal cells. It was demonstrated for the considered fully spin-polarized 
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Figure 3.7: The 3D projection of the angular part of the 10 S(sp 3 d 5 ) state Hartree-Fock 
node (with radial parts of orbitals identical for all spd orbitals). The projection has a 
topology of a single-sheet or double-sheet hyperboloid. The small cross shows the location 
of the nucleus while the spheres illustrate the electron positions. 

systems that the corresponding HF wave functions have the desired topology, i.e., two max- 
imal nodal cells. On the other hand, for partially spin-polarized and unpolarized systems 
the corresponding HF wave functions of the form ^ hf = det [ip a (r^jdetf^^r.,-)] lead to four 
nodal cells because there are two nodal cells for each independent spin subspace (assuming 
more than one electron in each spin subspace). To see how correlations eliminate a redun- 
dant nodal structure due to the artificial HF spin-up and spin-down separation we compare 
the nodal structure for HF and CI wave functions in two specific cases of spin-unpolarized 
states. 



3.5.1 HF and CI Nodes for Singlet State 1 S{p % 



The 3D projected nodal structure of 1 S(p®) unpolarized state is shown in Fig. 3.9 



We fixed positions of four electrons and sample for the nodal surface with a pair of electrons 
of opposite spins which are close to each other. The HF node consists of four nodal regions 
with a planar node for each spin subspace. The planar node is even exact for the case of 
three electrons in the same spin subspace, i.e., 4 S'(p 3 ). On the other hand, the CI wave 
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Figure 3.8: Projected Hartree-Fock node of 10 S(sp 3 d 5 ) of the majority spin valence elec- 
trons in Mn +2 ion. The Ne-core electrons are eliminated by pseudopotentials. Note the 
deformations from the radial parts of orbitals, including a small bubble detached from the 
rest of the surface (the right picture). For clarity, the positions of electrons have been 
omitted. 




Figure 3.9: The 3D projected HF and CI nodes of unpolarized 1 S(p e ) state. The HF node 
(left) consists of four nodal regions (two for each spin channel) while the CI node (right) 
exhibits only two nodal cells. Positions of four electrons (small spheres) are fixed, two 
different colors indicate opposite spins. The projected nodal surface is sampled with a pair 
of electrons of opposite spins which are positioned close together. 

function, which is close to the exact ground state of the given symmetry, leads to only two 
nodal cells. If we write *&ci = — e 2 ^HF + e^extit where ^ exc it includes single, double, 
and triple excitations, we find that e ~ 0.01. Introducing a small correlation via nonzero 
e makes a dramatic change to the topology of the HF nodal structure despite of a small 
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Figure 3.10: The 3D projected isosurfaces of HF and CI wave functions for 1 S(p e ) state. 
Different colors (dark and light) represent opposite signs. The HF wave function (left) has a 
discontinuity between regions of the same sign. In the CI wave function (right) the regions 
of the same sign are connected. A cut through (a view from top) reveals open channels 
(indicated by arrows) between regions of the same sign. 



total energy gain. Comparing the corresponding isosurfaces for HF and CI wave functions 



(Fig. 3.10) the structure of the wave functions seems to be very similar. However, a closer 
inspection reveals that regions of the same sign, which are disconnected in the HF wave 
function, are now connected by narrow channels in the CI wave function. 



3.5.2 HF and CI Nodes for N 2 Dimer State 

The projected nodal structure of the spin-unpolarized ground state 1 Y<~^ for N 2 



dimer is shown in Fig. 3.11 The He-core electrons are eliminated by pseudopotentials. The 
HF wave function with the separation of electrons into independent spin-up and spin-down 
subspaces forms fours nodal cells. In the CI case we can see that the HF nodes have been 
distorted and bended to build up channels, which connect regions of the same sign, i.e., two 
maximal nodal cells are obtained. 

In order to compare accuracy of the nodes we have carried out a fixed-node dif- 
fusion Monte Carlo simulations for N 2 dimer with core electrons eliminated by pseudopo- 
tentials pp. For the HF trial wave function (including the Jastrow factor, which does not 
change nodes) we have obtained the total energy -19.8395(7) H, which recovers 94% of cor- 
relation energy. For the CI trial function with more than 5000 determinants from double 
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Figure 3.11: The 3D projected nodes for N2 dimer. The HF nodes (left) are bended and 
distorted due to spin correlations into CI nodal surface (right) with just two nodal cells. 
Sampling of nodes is performed with a pair of electrons with opposite spins, which are close 
to each other; positions of other electrons (crosses) are fixed; small spheres indicate ions. 



excitations of 45 orbitals used and ccpV6Z basis (with up to / basis functions) the total 
energy reaches -19.870(5) H, which recovers 98% of correlation energy. The estimated exact 



energy is -19.8822 H. The Fig. 3.11 illustrates a good quality of the HF nodes except for 
the regions where two nodal surfaces cross each other. 



3.6 Conclusions 



We summarized the most recent knowledge about the general properties of nodes of 
fermionic wave functions. Further, we have investigated the nodes of atomic and molecular 
spin-polarized systems with one-particle states in s,p,d channels. We have studied cases 
with high symmetries, which enabled us to find exact nodes for several states with a few 
electrons (p 2 ,p 3 , 7r 2 ). Moreover, the projection of multi-dimensional manifolds into 3D space 
enabled us to study and characterize properties of nodes, in particular, their topologies for 
the Hartree-Fock wave functions. This analysis has provided useful insights and enabled us 
to formulate general transformation of one-particle coordinates using coordinate translation 
(backflow) and metric tensor to capture inhomogeneities and rotation symmetries. We test 
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a special form of a backflow coordinate transformation in Ch. [5] of the thesis. Finally, 
we have demonstrated that the more accurate CI wave functions for two specific cases of 
spin-unpolarized states have the nodal structure of only two maximal nodal cells. 



56 



Chapter 4 



Pfaffian Pairing Wave Functions 

Sections of this chapter also appeared in: 
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4.1 Introduction 

The key challenge for successful application of fixed-node DMC is to develop meth- 
ods, which can eliminate the fixed-node bias or at least make it smaller than experimental 
error bars for the given quantity. This is a difficult task, once we realize that the fermion 
nodes, which are subject of Ch. |3j are complicated high-dimensional manifolds determined 
by the many-body effects. So far, improvement in the accuracy of trial wave functions has 
proved to be one realistic approach to finding better approximations for the nodes. This 
approach has an additional benefit in forcing us to think about the relevant correlation 
effects and their compact and computationally efficient description. 

The commonly used QMC trial wave functions have the Slater-Jastrow form, which 
can be written as = ^A^ x p[U C orr], where a is the antisymmetric part while U corr de- 



scribes the electron-electron and higher-order correlations as described in Sec. 2.4.2 The 
antisymmetric component is typically one or a linear combination of several Slater deter- 
minants of one-particle orbitals such as a configuration interaction expansion introduced in 



Sec. 1.3. To overcome the limit of one-particle orbitals the two-particle or pair orbital has 
been suggested. In condensed systems one such example is the Bardeen-Cooper-Schrieffer 
(BCS) wave function, which is an antisymmetrized product of singlet pairs. The singlet 
pair is sometimes referred to as geminal and the resulting wave function as the antisym- 
metrized geminal product (AGP). It has been recently used to calculate several atoms and 
molecules as well as superfluid Fermi gases |94[ [95) 1105] . The results show promising gains 
when compared to the single-determinant Hartree-Fock (HF) wave functions, nevertheless, 
in partially spin-polarized systems the improvements are less pronounced due to the lack of 
pair correlations in the spin-polarized subspace [94, 95j. The spin-polarized (triplet) pairing 
wave functions lead to Pfaffians (instead of determinants) and have been mentioned a few 
times before and applied to model systems |106[ I107[ 1108] . 



In this chapter, we further develop the idea from Sec. 2.4.2 in which we have 



proposed the description of electronic systems by a novel generalized pairing wave function 



in the Pfaffian form. This chapter is organized as follows: In Sec. 4.2 we present a set of 



the key mathematical identities and formulas for Pfaffians, some of them derived for the 



first time. In Sec. 4.3, we establish the connection of a generalized Pfaffian pairing wave 
function to BCS and HF wave functions. The resulting Pfaffian wave functions are tested on 
atomic and molecular systems in variational and fixed-node diffusion Monte Carlo methods 
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as described in Sec. 4.4.1 In Sees. 4.4.2 and 4.4.3 we investigate also a generalizations 
to linear combinations of Pfafflans and to antisymmetrized independent singlet pairs and 
compare the results from the point of view of recovered energies and compactness of the 



wave functions. Finally, in Sec 4.4.4 we analyze the fermion nodes for some of the wave 
functions and point out the topological differences between HF, Pfaffian and an essentially 
exact wave functions for a given test example. 



4.2 Algebra of Pfaffians 



4.2.1 Definitions 

First introduced by Arthur Cayley in 1852 |109| . the Pfaffian is named after 
German mathematician Johann Friedrich Pfaff. Given a 2n x 2n skew-symmetric matrix 
A = [a>ij] , the Pfaffian of A is defined as antisymmetrized product 



pf [A] = A[a lt2 a 3A . . . a 2n -i,2n] 



(4.1) 



where the sum runs over all possible (2n— 1)!! pair partitions a = (i 2 , J2), • • • , (in, jn))} 

of {1, 2, ... , 2n} with < j^. The sign of permutation associated with the partition a is 
denoted as sgn(a). The Pfaffian for a matrix of odd order equals to zero. The following 
example gives Pfaffian of a A(4 x 4) skew-symmetric matrix 








a\2 


ai3 




pf 


-«12 





023 


024 


-Ol3 


-023 





034 




-au 


-024 


-034 






012034 — 013024 + ai4<223- 



(4.2) 



It can be also evaluated recursively as 

2n 



P f [^] = Yl ttl 'J Yl S S n ( a lj) a h,h a i 2 j2 ■ ■ ■ a in-l,3n-l 

2n 

i=2 



(4.3) 
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where a\ j is partition with ik,jk 1>J an d Pc{ a ij) is defined as Pfaffian cofactor of aij. 
The cofactor for an element aj k is given by a formula, 

P c (o i)fc ) = (-iy +k+1 pf [A(j, k; j, k)} (4.4) 

where the matrix A(j, k;j, k) has the rank 2(n — 1) x 2(n — 1) and is obtained from A by 
eliminating j and k rows and columns. 



4.2.2 Calculation of a Pfaffian 



There exist several identities involving Pfaffians and determinants. For any 2n x 2n 
skew-symmetric matrix A and arbitrary matrices B(2n x 2n) and M{n x n) we have the 
following relations: 



pf 



P fL4 T ] = (-1)VL4] 
pi[A} 2 = det[A] 

Ai 
A 2 

pf [BAB' 1 '] = det[B]pf[A] 



pf[A!]pf[A 2 ] 



pf 



Key ideas of respective proofs: 





-M 1 



M 



n(n — l) 

(-lJ-V^detlAf] 



(4.5a) 
(4.5b) 

(4.5c) 

(4.5d) 

(4.5e) 



(4.5a) Each permutation contains product of n pairs resulting in an overall (— l) n factor. 



(4.5b) This is a well-known Cayley's relationship between the Pfaffian and the determinant 
of a skew-symmetric matrix. Since it has been proved many times before in variety 
of ways |110| 1112] we do not give this proof here. Using this relation we rather 
prove a more general version of Cayley's identity |112] in App. [B] which we were not 
able to find anywhere else except in the original Cayley's paper |112j . 



(4.5c) Use the expansion by Pfaffian cofactors. 



(4.5d) By squaring (4d), using Eq. (4.5b), and taking the square root one finds pi[BAB T ] = 
±det[-B]pf [A]. Substituting the identity matrix / for B one finds + to be the correct 



sign. 
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(4.5e) Assume 



B 



M 
/ 



and A 



I 
-I 



in Eq. (4.5d). The overall sign is given by value of pf[.A]. 



The identities listed above imply several important properties. First, Eqs. (4.5d| 



and (4.5e) show that every determinant can be written as a Pfaffian, but on the contrary, 



only the absolute value of Pfaffian can be given by determinant [Eq. (4.5b)]. The Pfaffian 
is therefore a generalized form of the determinant. Second, by substituting suitable matri- 



ces [1 13j for M in Eq. (4.5d) one can verify the following three properties of Pfaffians |114j . 
similar to the well-known properties of determinant. 

(a) Multiplication of a row and a column by a constant is equivalent to multiplication of 
Pfaffian by the same constant. 

(b) Simultaneous interchange of two different rows and corresponding columns changes 
the sign of Pfaffian. 

(c) A multiple of a row and corresponding column added to to another row and corre- 
sponding column does not change the value of Pfaffian. 

It is also clear that any skew-symmetric matrix can be brought to block-diagonal form by 



an orthogonal transformation. Recursive evaluation [Eq. (4.3)] implies that the Pfaffian of 
block-diagonal matrix is directly given by 



pf 





-Ai 



Ai 






-A 2 



A 2 






-A, 



An 




A1A2 . . . A n . 



(4.6) 



Therefore by employing a simple Gaussian elimination technique with row pivoting (see 
App. |C| we can transform any skew-symmetric matrix into block-diagonal form and obtain 
its Pfaffian value in 0(n 3 ) time. 
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However, in QMC applications, one often needs to evaluate the wave function after 
a single electron update. Since Cayley |112] showed (for proof see App. pi) that 



det 
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■ 0-2n 
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we can relate the Pfaffian of original matrix pf [A] to the Pfafnan of a matrix with updated 
first row and column pf [B] using the inverse matrix A -1 in only 0(n) operations by 

det[A] EjhjAj, 1 



pf[B] 



pf[A] 



pflA^hjAj*. 



(4i 



The second part of Eq. (4.8) was obtained by taking advantage of the identity in Eq. (4.5b). 
Slightly more complicated relation between pf [A] and pf [B] can be derived if one considers 
simultaneous change of two separate rows and columns, which represents the two electron 
update of a wave function. 



4.2.3 Gradient and Hessian of Pfafnan 



In the case of linear dependence of matrix elements Aona set of parameters {c}, 
one can derive the following useful relations: 



1 dpf[A] 1 



and 



1 d 2 pf[A] _1 
pf[A] dadcj ~4 



tr 



. x dA 
da 



tr 



' A -iOA' 




A -iOA- 
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(4.9) 



(4.10) 
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- -tr 
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- A -idA A _ x dA 
da da 



where A is again the inverse of A. 
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4.3 Pairing Wave Functions 

In order to contrast the properties of pair wave functions with the ones build from 
one-particle orbitals we will first recall the well-known fact from the Hartree-Fock theory. 
The simplest antisymmetric wave function for N electrons constructed from one-particle 
orbitals is the Slater determinant 

Vhf = det[£ fe (ri,CTi)] = det[<p k (i)]; i,k = l,...,N, (4.11) 

where tilde means that the one-particle states depend on both space and spin variables. 
Clearly, for N electrons this requires N linearly independent spin-orbitals, which form an 
orthogonal set in canonical HF formulation. 

Let us now consider the generalization of the one-particle orbital to a two-particle 
(or pair) orbital 4>(i,j), where tilde again denotes dependence on both spatial and spin 
variables. The simplest antisymmetric wave function for 2N electrons constructed from the 
pair orbital is a Pfaffian 

* = A[j>(l, 2), 0(3,4) . . . 4>(2N - 1,2N)\ = pf[4>(i,j)}. (4.12) 



The antisymmetry is guaranteed by the definition (4.1), since the signs of pair partitions 
alternate depending on the parity of the corresponding permutation. The important differ- 
ence from Slater determinant is that in the simplest case only one pair orbital is necessary. 
(This can be generalized, of course, as will be shown later.) If we further restrict our de- 
scription to systems with collinear spins, the pair orbital 0(rj, af, Tj, o~j) for two electrons 
in positions rj and rj and with spins projections crj and <7j can be expressed as 

faiMTj,*-) = <Ki,j)(<Ti<T3\[\ TI) " | IT)]/V2 (4.13) 

+ x n (hj)(<Ti<Tj\ TT) 

+ x n (i,j)M[\ U) + UT)]/V2 

+ X U (i,j){<nc7j\ ID- 
Here = </>(rj,r,-) is even while x > X an d X are °dd functions of spatial coordi- 



nates. In the rest of this section we will discuss special cases of the wave function (4.12). 
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4.3.1 Singlet Pairing Wave Function 

Let us consider the first 1, 2, iV electrons to be spin- up and the rest N + l, ...,2N 
electrons to be spin-down and allow only 0(rj,rj) in </>(rj, o~f, Tj, Gj) to be non-zero. Using 



the Pfaffian identity [Eq. (4.5e)] we can write the wave function for N singlet pairs, also 



known as the BCS wave function (or AGP), in the following form 



^BCS = Pf 



* n 

_$UT 



det[* n ], (4.14) 



which is simply a determinant of the N x N matrix = [4>(i,j)] as was shown previ- 
ously [TT5l[TT6] . 

It is straightforward to show that the BCS wave function contains the restricted 
HF wave function as a special case. Let us define the Slater matrix C = \fi{j)\ where {tpi} 
is a set of HF occupied orbitals. Then we can write 

*hf = det[C]det[C] = det[CC T ] = det[*{J>], (4.15) 

where 

N 

(*h>)vj = <i>HF(hj) = Y J ( Pk(i)Vk(j)- (4.16) 

fc=l 

On the other hand, we can think of the BCS wave function as a natural generalization of 
the HF one. To do so we write the singlet pair orbital as 

M>N 

4>{i,j) = S km {i)m{j) = <p(i)S<p(j), (4.17) 

k,l 

where the sum runs over all M (occupied and virtual) single-particle orbitals and S is some 
symmetric matrix. Therefore, we can define one-particle orbitals, which diagonalize this 
matrix and call them natural orbitals of a singlet pair. 

The BCS wave function is efficient for describing systems with single-band corre- 
lations such as Cooper pairs in conventional BCS superconductors, where pairs form from 
one-particle states close to the Fermi level. 

4.3.2 Triplet Pairing Wave Function 

Let us assume, in our system of 2N electrons, that the first Mi electrons are spin- 
up and remaining M<i = 2N — Mi electrons are spin-down. Further, we restrict M\ and M2 
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to be even numbers. Then by allowing only x^(i, j) and X m (4.13) to be non-zero, 



we obtain from expression (4.12) by the use of Eq. (4.5c 



111 



(4.18) 





£ u 

where we have introduced M\ x M\{M.2 x M2) matrices ^TT(J-l) = [x^^^(i, j)] . To our 
knowledge, this result was never explicitly stated and only the weaker statement that the 
square of wave function simplifies to a product of determinants has been given |115[ 1116] . 

The connection to a restricted HF wave function for the above state can be again 
established as follows. In accord with what we defined above, det[(C)^^] are spin-up(-down) 



Slater determinants of some HF orbitals {<fi}. Then, by taking advantage of Eq. (4.5e) we 
can write 



^hf = det[C T ]det[C^] 



(4.19) 



pf[CUiCt T ]pf[CU 2 C^ J 



pfL^pfL/y 

given Ai and A2 are some skew-symmetric non-singular matrices. In the simplest case, 
when A\ and A2 have block-diagonal form ( |4.6[ ) with all values Aj = 1, one gets 



The pair orbitals can be then expressed as 



h3 



M 1 (M 2 )/2 

E ^ 

k=l 



(4.20) 



(4.21) 



2fc-l 



(i)<P2k(j) - <P2k-l{j)<P2k(i))- 



Similarly to the singlet pairing case, one can also think of the triplet pairing wave function 
as a natural generalization of the HF one. To do so we write the triplet pair orbitals as 



,TT(U) 



M>M 1 (M 2 ) 



(i)A TT(U)^ (j)) 



(4.22) 



where again the sum runs over all M (occupied and virtual) single-particle orbitals and 
A.TT(H) are some skew-symmetric matrices. Therefore, we can define one-particle orbitals, 
which block-diagonalize these matrices and call them natural orbitals of a triplet spin-up- 
up( down- down) pair. 
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4.3.3 Generalized Pairing Wave Function 

Let us now consider a partially spin-polarized system with unpaired electrons. In 



order to introduce both types of pairing we allow x (hj)> X (*>j) an d 4>(hj) i n (4.13) 
to be non-zero. However, we omit the X^ihj) term. Then our usual ordered choice of 
electrons labels with all spin-up electrons first and remaining electrons spin-down enables 



us to directly write from (4.12) the singlet-triplet-unpaired (STU) orbital Pfaffian wave 
function [14J 



STU 



Pf 



-<f>1lT fill 



ip I 



AT 



(4.23) 



where the bold symbols are block matrices or vectors of corresponding orbitals as defined 



in Sections 4.3.1 and 4.3.2 and T denotes transposition. For a spin-restricted STU wave 
function the pair and one-particle orbitals of spin-up and -down channels would be identical. 

The Pfaffian form can accommodate both singlet and triplet pairs as well as one- 
particle unpaired orbitals into a single, compact wave function. The correspondence of STU 
Pfaffian wave function to HF wave function can be established in a similar way to the pure 
singlet and triplet pairing cases. 



4.4 Results 



The Pfaffian wave functions were used in QMC calculations by variational and 
fixed-node diffusion Monte Carlo. As we have mentioned earlier, the VMC trial wave 
function is a product of an antisymmetric part $?a times the Jastrow correlation factor 



*y M c(R) = *A(R)exp[J7 corr ({r i:7 }, {ru}, {rji})}, 



(4.24) 



where U corr depends on electron-electron, electron-ion and, possibly, on electron-electron- 



ion combinations of distances as described in Sec. 2.4.2 For the antisymmetric part we have 
used ^ a = ^ hf and \t a = ^ stu- Some tests were also done with ^ a = ^ bcs to compare 
with recent results |94| [§5] . The pair orbitals were expanded in products of a one-particle 



orbital basis [94 according to Eqs. (4.17) and (4.22). The expansions include both occu- 



pied and unoccupied (virtual) one-particle orbitals. The one-particle atomic and molecular 
orbitals used were either Hartree-Fock orbitals or natural orbitals [39] from CI correlated 
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calculations. Typically, we used about 10 virtual orbitals. The natural orbitals produced 
better and more systematic results than the HF ones. The pair orbital expansion coeffi- 
cients were then optimized in VMC by minimizations of energy, variance or a combination 



of energy and variance (for details, see Sec. 2.6.4). The optimization procedure requires the 



calculation of gradient and the Hessian of the wave function according to Eqs. (4.9) and 



(4.10). 



4.4.1 Single Pfaffian Calculations 

We have applied these developments to several first row atoms and dimers (see 



Table 4.1 and Fig. 4.1 ). We used pseudopotentials to eliminate the atomic cores [HE], while 
the previous all-electron calculations with BCS wave functions |94| [95j produced percentages 
of the correlation energies in accordance with our BCS wave functions calculations. 

Perhaps the most striking result is a systematic percentage of recovered correlation 



energy on the level of 94-97% in DMC method for the STU wave functions (see Table 4.1 



and Fig. 4.1 ). Another significant result is that in general the triplet contribution for these 
single Pfaffian STU wave functions are small, with the only exception being nitrogen atom, 
where we see a gain of additional 1% in correlation energy when compared to a trial wave 
function without triplet pairs. We believe, this is due to the fact, that ground state of 
nitrogen atom has a quartet spin state and therefore the highest spin polarization from all 
studied cases. However, only future tests with the presence of non-zero X (hj) m (4.13) 



will determine the full extent of the triplet contribution to the correlation effects. Given 
the pair orbitals have been optimized using VMC method, it is natural that the relative 
gains in correlation energy with respect to the HF wave functions are larger on the level of 
VMC calculations than in DMC calculations. Overall, the single Pfaffian form is capable 
of capturing near-degeneracies and mixing of excited states for both spin-polarized and 
unpolarized systems. 

4.4.2 Multi-Pfaffian Calculations 

To test the limits of the Pfaffian functional form, we have proposed a simple 
extension: the multi-Pfaffian (MPF) wave function of the form 

^mpf = wipt [x[\ x[ l , <£i, Vl] + w 2 pf [xl\ X l 2 l , 02, y%] + ■ ■ ■ , (4.25) 
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Figure 4.1: Correlation energies obtained by QMC methods with the different trial wave 
functions: VMC and fixed-node DMC with HF nodes (HF) and STU Pfaffian nodes (PF). 
The lower plot shows the fixed-node DMC correlation energy gains over HF nodes for BCS 
and STU Pfaffian wave functions. The statistical error bars are of the symbol sizes or 
smaller. Except for the Be atom all the calculations used the same pseudopotentials [TJ [2]. 



where Wi denotes the weight of ith Pfaffian. In order to improve upon the wave function 



with single STU Pfaffian, the additional terms in wave function (4.25) have to contain 
some new excitations not previously present. As an example of this form, we apply it to the 



carbon pseudo-atom. The pair orbitals, Eqs. (4.17) and (4.22 ), for this system are expanded 
in the basis of HF occupied orbitals 2s, 2p x and 2p y . The choice of singlet </>i(l,2) = 
2s(l)2s(2) = 0i[2s,2s] and spin-up spin-up triplet xI T (l, 2) = 2p x (l)2p y (2)-2p y (l)2p x (2) = 
x\^[2p x , 2p y ] pair orbitals (the other functions are taken to be zero) then gives pf[xi j 4>i] = 
^hf[2s^, 2p x , 2p y \. However, one can construct the equivalent combinations of pairs as 
4>2[2s,2p x ], X2^[2s, 2p y \ and 4>3[2s, 2p y ], X3^[2s,2p x ]. We can therefore include all three 
Pfaffians into our <!/ mpf and further optimize all the pairing functions in VMC method on 
the space of occupied and virtual orbitals. 

Since each pair orbital in ^ mpf contains on the order of M 2 pairing coefficients, 
M being the total number of one particle orbitals involved, we limit our expansions to only 
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few Pfaffians. However, this can be improved by factor M, if we diagonalize the coefficient 
matrices. In practice, given the optimization routine in VMC method can safely minimize 
on the order of a few hundred coefficients at the same time, we end up doing several 
partial optimizations. To minimally disrupt the weights in partially expanded MPF wave 
function, any new Pfaffians are added in pairs, each initially set to the HF wave function 
with opposite sign for zero net contribution. Then the Pfaffian pair is re-optimized on the 
set of all single-particle orbitals. Besides already optimized STU Pfaffian we have an even 
number of additional Pfaffians, which explains the overall odd number of Pfaffians in our 
MPF expansions. 



The results in Table 4.2 show that for the atomic systems our MPF wave functions 
are able to recover close to 99% of correlation energy. Furthermore, comparison with the 
CI results demonstrates it is possible to obtain similar quality wave functions with corre- 
sponding improvements of the fermion nodes at much smaller calculational cost. However, 



for the diatomic cases (see Table 4.3), only very limited gain over single STU Pfaffian wave 
function correlation energies were achieved for MPF wave functions with few Pfaffians. We 
therefore conclude that for obtaining significantly larger gains in correlation the molecular 
wave functions require much larger expansions. 



4.4.3 Antisymmetric Independent Pairs Wave Function 

We have also tested the fully antisymmetric independent pairs (AIP) wave function 
which introduces one pair orbital per each electron pair. For system of 2N fermions in singlet 
state the AIP wave function can be written as 

q, AIP = A[<k(l, 2), 02(3, 4), . . . , 4> N (2N - 1, 2N)) (4.26) 

= ^pf[4(l,2),4(3,4),...,0i JV (2iV-l,2iV)], 
p 

where the last equation corresponds to the sum over all N\ possible permutations of N 
different pair orbitals <pi for each Pfaffian. This wave function is closely related to the wave 
function of an antisymmetrized product of strongly orthogonal geminals |118j . The results 
for C2 and N2 dimers using AIP wave function are given in Table |4~3} 

Consideration of independent pairs results in an exponential increase of a number 
of Pfaffians. However, captured correlation energy is on the level of small MPF expansion, 
and significantly less than CI with re-optimized weights using the same one-particle orbitals. 
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Figure 4.2: A three-dimensional cut through the fermion node hypersurface of oxygen atom 
obtained by scanning the wave function with a spin-up and -down (singlet) pair of electrons 
at the equal positions, while keeping the rest of electrons at a given VMC snapshot positions 
(small green spheres). Nucleus is depicted in the center of the cube by the blue sphere. The 
three colors (from left to right) show nodes of: Hartree-Fock (red/dark gray); Multi-Pfaffian 
nodes (orange/medium gray); and the nodes of the CI wave function (yellow/light gray) in 
two different views (upper and lower rows). The CI nodal surface is very close to the exact 
one (see text). The HF node clearly divides the space into four nodal cells while Pfaffian 
and CI wave functions partitioning leads to the minimal number of two nodal cells. The 
changes in the nodal topology occur on the appreciable spatial scale of the order of 1 a.u. 

This suggests that to achieve more correlation energy in larger systems we have to go beyond 
double pairing. 

4.4.4 Nodal Properties 

As we have already mentioned in Ch. [3] the fermion node manifold is defined by an 
implicit equation ^(i?) = and for N electrons it is a (3iV — l)-dimensional hypersurface. 
With exception of a few exact cases, the nodes of trial wave functions introduce bias into 
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the fixed-node DMC energies. Recently a number of authors have reported improvement in 
nodal structure of trial wave functions [Til [23l 1241 l92l 193] [Ml 195] . 

The effect of pairing correlations on nodes can be highlighted by direct comparison. 



The Fig. 4.2 shows the example of nodal structure of oxygen atom. Here we compare the 
nodal surfaces of HF (no pairing), MPF Pfaffian (STU pairing) and a high accuracy CI 
wave function with more than 3000 determinants, which gives essentially exact fermion 
nodes (i.e., 99.8(3)% of correlation energy in fixed-node DMC). 

It is clear that the changes in the nodal surfaces are significant, the most important 
one being the elimination of artificial four nodal cells resulting from the independence of 
spin-up and -down channels in HF. The Pfaffian smooths-out the crossings and fuses the 
compartments of the same sign into the single ones. These topology changes therefore 
lead to the minimal number of two nodal cells, an effect observed in correlated context 
previously |83[ 1851 1861 l9T| [96] . However, the nodes of the Pfaffian wave functions could be 
further improved if the scheme for direct optimization of nodes of trial wave functions were 
used |105[ I119j . Another interesting result from our work is that despite such a substantial 
change in the nodal structure the amount of missing correlation energy is still non-negligible. 



4.5 Conclusions 

To summarize, we have proposed Pfaffians with singlet pair, triplet pair and un- 
paired orbitals as variationally rich and compact wave functions. They offer significant and 
systematic improvements over commonly used Slater determinant-based wave functions. 
We have included a set of key mathematical identities with proofs, which are needed for the 
evaluation and update of the Pfaffians. We have also shown connections of HF and BCS 
(AGP) wave functions to more general Pfaffian wave function. We have further demon- 
strated that Pfaffian pairing wave functions are able to capture a large fraction of missing 
correlation energy with consistent treatment of both spin-polarized and unpolarized pairs. 
We have explored multi-Pfaffian wave functions which enabled us to capture additional 
correlation. While for atomic systems the results are comparable to large-scale CI wave 
functions, molecular systems most probably require much larger multi-Pfaffian expansions 
than we have explored. As another test of the variational potential of pairing we have em- 
ployed the fully-antisymmetrized independent pairs wave function in Pfaffian form and we 
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have found that it does not lead to additional gains in correlation energy. We therefore con- 
clude that more general functional forms together with more robust large-scale optimization 
methods might be necessary in order to obtain further improvements. The gains in correla- 
tion energy for Pfaffians come from improved fermion nodes, which are significantly closer 
to the exact ones than the HF nodes, and exhibit the correct topology with the minimal 
number of two nodal cells. 
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Table 4.1: Total energies for C, N and O atoms and their dimers with amounts of the 
correlation energy recovered in VMC and DMC methods with wave functions as discussed 
in the text. Unless noted otherwise, the numbers in parentheses are the statistical errors in 
the last digit from corresponding QMC calculation. Energies are in Hartree atomic units. 
For C, N, O atoms we used the correlation energies by Dolg [4J(0.1031, 0.1303, 0.1937 H). 
For the estimation of correlation energies of dimers we needed accurate HF energies at 
experimental distances [5] and the estimated exact total energies. Each exact total energy 
was estimated as a sum of total energies of constituent atoms minus experimental binding 
energy [5j El |7] adjusted for experimental zero-point energy [7] . 



Method/ WF 


C 


Ecorr [%] 


N 







E C orr[%] 


HF 


-5.31471 





-9.62892 





-15.65851 





VMC/HF 


-5.3939(4) 


76.8(4) 


-9.7375(1) 


83.3(1) 


-15.8210(6) 


83.9(3) 


VMC/BCS 


-5.4061(2) 


88.6(2) 


-9.7427(3) 


87.3(2) 


-15.8250(3) 


86.0(2) 


VMC/STU 


-5.4068(2) 


89.3(2) 


-9.7433(1) 


87.8(1) 


-15.8255(3) 


86.2(2) 


DMC/HF 


-5.4061(3) 


88.6(2) 


-9.7496(2) 


92.6(2) 


-15.8421(2) 


94.8(1) 


DMC/BCS 


-5.4140(2) 


96.3(2) 


-9.7536(2) 


95.7(2) 


-15.8439(4) 


95.7(2) 


DMC/STU 


-5.4139(2) 


96.2(2) 


-9.7551(2) 


96.8(1) 


-15.8433(3) 


95.4(2) 


Est. /Exact 


-5.417806 


100 


-9.759215 


100 


-15.85216 


100 


Method/ WF 


c 2 


E corT [%] 


N 2 


E corr [%] 


o 2 


E corr [%] 


HF 


-10.6604 





-19.4504 





-31.3580 





VMC/HF 


-10.9579(4) 


72.9(1) 


-19.7958(5) 


80.0(1) 


-31.7858(6) 


79.6(1) 


VMC/BCS 


-11.0059(4) 


84.7(1) 


-19.8179(6) 


85.0(1) 


-31.8237(4) 


86.7(1) 


VMC/STU 


-11.0062(3) 


84.8(1) 


-19.821(1) 


85.8(2) 


-31.8234(4) 


86.6(1) 


DMC/HF 


-11.0153(4) 


87.0(1) 


-19.8521(3) 


93.0(1) 


-31.8649(5) 


94.3(1) 


DMC/BCS 


-11.0416(3) 


93.5(1) 


-19.8605(6) 


94.9(1) 


-31.8664(5) 


94.6(1) 


DMC/STU 


-11.0421(5) 


93.6(1) 


-19.8607(4) 


95.0(1) 


-31.8654(5) 


94.4(1) 


Est./Exact c 


-11.068(5) a 


100.0(10) 


-19.8825(6) b 


100.0(1) 


-31.8954(l) b 


100.0(1) 



a There is rather large discrepancy in the experimental values of C2 binding energy (141.8(9) 143(3) [7] 
and 145.2(5) kcal/mol 6J. For the estimation of exact energy we have taken the average value of 143(3) 
kcal/mol. 

b Experimental binding energies taken from ref. [5]. 

c The error bars on estimated exact total energies are due to experiment. 
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Table 4.2: Percentages of correlation energies recovered for C, N and O atoms by VMC and 
DMC methods with wave functions as discussed in the text. The corresponding number of 
Pfaffians or determinants n for each wave function is also shown. For details, see caption 
of Table 



4.1 



Method/ WF 


n 


C 


n 


N 


n 





VMC/MPF 


3 


92.3(1) 


5 


90.6(1) 


11 


92.6(3) 


VMC/CP 


98 


89.7(4) 


85 


91.9(2) 


136 


89.7(4) 


DMC/MPF 


3 


98.9(2) 


5 


98.4(1) 


11 


97.2(1) 


DMC/CP 


98 


99.3(3) 


85 


98.9(2) 


136 


98.4(2) 



The determinantal weights were taken directly from CI calculation without re-optimization in VMC. 



Table 4.3: Total energies for C2 and N2 dimers with amounts of correlation energy recovered 
in VMC and DMC methods with wave functions as discussed in the text. Energies are in 
Hartree atomic units. The corresponding number of Pfaffians or determinants n for each 
wave function is also shown. For details, see caption of Table 



4.1 



Method/ WF 


11 


c 2 


E corr [%] 


11 


N 2 


E corr [%] 


VMC/MPF 


5 


-11.0187(2) 


87.8(1) 


5 


-19.8357(3) 


89.2(1) 


VMC/AIP 


4! 


-11.0205(4) 


88.3(1) 


5! 


-19.8350(3) 


89.0(1) 


VMC/CP 


148 


-11.0427(1) 


93.7(1) 


143 


-19.8463(9) 


91.6(2) 


DMC/MPF 


5 


-11.0437(4) 


94.0(1) 


5 


-19.8623(5) 


95.3(1) 


DMC/AIP 


4! 


-11.0435(7) 


94.0(2) 


5! 


-19.8611(3) 


95.0(1) 


DMC/CI a 


148 


-11.0573(2) b 


97.3(1) 


143 


-19.875(2) 


98.3(5) 



The determinantal weights were re-optimized in the VMC method. 

Recently, Umrigar et.al. [93] published very accurate DMC result for fully optimized CI wave function 
with up to 500 determinants for C2 molecule. The resulting well-depth of his calculation is 6.33(1) 
eV, which is only 0.03 eV form estimated exact value of Ref. |117| . The well-depth resulting from our 
DMC/CI energy of -11.0573(2) H equals to 6.03(1) eV. 
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Chapter 5 



Backflow Correlations in Slater 
and Pfaffian Wave Functions 

5.1 Introduction 

Another approach for improvement of the nodal accuracy of variational trial wave 
functions is to employ backflow correlations first introduced by Feynman and Cohen [16] for 
liquid 4 He. Since then, a number of authors [EJ [EJ 122 EH [22] showed that the backflow 
correlations are helpful also in fermionic systems. Recently, the backflow transformation 
of electron coordinates has been applied to chemical (inhomogeneous) systems |23| I24j . In 
this chapter, we report on application of backflow transformation, already introduced in 
the Sec. |2.4.2| to two simple but nontrivial testing cases: carbon atom and its dimer. It is 
the first application of backflow correlations to the multi-determinantal and Pfaffian wave 
functions. 
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5.2 Backflow Wave Function Form 



As already discussed in Sec. 2.4.2, our trail wave function has the form 



* r (R) = tf A (X) x exp[U corr (R)}, (5.1) 

where X = (xi,...,xjv) represents some quasi-particle coordinates dependent on all N 
electron positions R. Further, $ a is either (multi)-determinantal or Pfaffian wave function 
and U corr is the Jastrow correlation factor both defined in the previous chapters. 
The quasi-coordinate of ith electron at position rj is given as 

Xj = Ti + £j(R) 

= r, + £»(R) + £f (R) + r n (R)> (5-2) 

where is the ith electron's backflow displacement divided to the contributions from one- 
body (electron-nucleus), two-body (electron-electron) and three-body (electron-electron- 
nucleus) terms. They can be further expressed as 

^ n (R) = J2x(ni)ru (5.3) 

*f (R) = J>(r y )r« ( 5 - 4 ) 

£f eri (R) = ^^[^(r^,^/,^/)^ + w 2 (ry,r i /,r J -/)ri/], (5.5) 

where rjj = rj — r,- and rj/ = rj — 17. The Xj u an d w\ with ^2 terms are similar to one, 
two and three-body Jastrow terms and are further expanded as 

X(r) = ^2c k a k (r), (5.6) 

k 

u(r) = Y / dkb k (r), (5.7) 

k 

wi,2(nj,ru,rji) = ^2gklmak(ni)ai(rjj)b m (nj). (5.8) 

klm 

The one dimensional basis functions {0} and {&} are chosen as Gaussians with the center 



in origin to preserve the electron-electron [Eqs. (2.26) and (2.27)] and electron-nucleus 



[Eq. (2.25 )] cusp conditions. The set of variational parameters {c}, {d} and {g} is minimized 



with respect to mixture of energy and variance in NEWTON_OPT method from Sec. 2.6.4 
In addition, all electron-electron coefficients ({d k } and {gkim} with fixed k and I) are allowed 
to be different for spin-like and for spin- unlike electron pairs. 
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PP carbon atom 



PP carbon dimer 
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Figure 5.1: Slater- Jastrow (SJ), Pfaffian-Jastrow (PF) and Cl-Jastrow (CI) wave functions 
with backflow (BF) correlations for carbon pseudopotential (He core) atom (left) and dimer 
(right) tested in VMC and DMC methods. Upper figure: Percentages of correlation energy 
versus a number of backflow parameters. Notation: 2B for all electron-nucleus and electron- 
electron terms, 3B for all electron-electron-nucleus terms only and 23B for all terms together. 
Lower figure: Variance of the local energy versus a number of terms. 

5.3 Results 

We test the above backflow correlation function on carbon pseudopotential (He 



core) atom and its dimer. The main results are summarized in Fig. 5.1 and can also be 



found in full detail in Tables 5.1 and 5.2 The backflow correlations are able to capture 
additional few percent of correlation energy for both Slater-Jastrow and Pfaffian-Jastrow 
wave functions. Another important feature of backflow is 20-30% decrase in variances of 
local energy with respect to the wave functions without backflow correlations. The gains are 
systematic with increasing number of parameters, however we do not find the three-body 
terms as important as previous study |24| . This can be attributed to the fact that their 
main testing case was an all-electron system where the three-body correlations are known 
to be more significant. 

In order to gain further insight into the action of backflow transformation on 
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012345678 012345678 
r [a.u.] r [a.u.] 



Figure 5.2: Optimized electron- nucleus x{ r ) an d electron-electron u{r) two-body backflow 
functions of carbon PP atom (left) and dimer (right) without the presence of three body 
term versus distance r. The functions shown are for Slater-Jastrow (SJ) and Pfaffian- 
Jastrow (PF) wave functions with backflow (BF) correlations. Note that = u^, i.e., 
this function is the same for all spin-like electron-electron pairs and similarly u'^ = is 
the same for all spin-unlike electron-electron pairs. 



electron coordinates we plot the electron-nucleus and electron-electron two-body backflow 
functions optimized in VMC method without the presence of three-body terms. As it is 



immediately visible from Fig. 5.2 the spin-unlike electron-electron functions are almost 
order of magnitude larger than electron-nucleus and spin-like electron-electron functions. 
They are also characterized with well defined shape, which is almost independent of type 
of trial wave function. 



5.4 Conclusions 

We have presented the first application of Pfaffian and multi-determinantal wave 
functions with backflow correlations to chemical systems. Results for two testing cases of 
carbon pseudo atom and its dimer show promising gains in correlations energies, decreases 
in variances and improvements in the nodal structure. However, it all comes at the addi- 
tional computational cost of calculation of the backflow displacement and the simultaneous 
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Table 5.1: Slater- Jastrow (SJ), Pfaffian-Jastrow (PF) wave functions with backflow (BF) 
correlations for carbon pseudo atom are tested in VMC and DMC methods. Notation is 
the same as in Fig. 



5.1 



Method WF 




N u 


N W1 


N W2 




E[H] 


a 2 [H 2 ] 


E corr .[%] 


HF S 












-5.31471 




0.0 


VMC SJ 












-5.3990(1) 


0.0677 


81.8(1) 


SJBF2B 


11 


22 






33 


-5.4011(2) 


0.0544 


83.8(2) 


SJBF3B 






128 


128 


256 


-5.4023(3) 


0.0504 


85.0(3) 


SJBF23B 


4 


8 


128 


128 


268 


-5.4020(2) 


0.0498 


84.7(2) 


PF 












-5.4083(2) 


0.0626 


90.8(2) 


PFBF2B 


11 


22 






33 


-5.4097(1) 


0.0427 


92.1(1) 


PFBF23B 


4 


8 


128 


128 


268 


-5.4107(1) 


0.0411 


93.1(1) 


DMC SJ 












-5.4065(3) 




89.0(3) 


SJBF2B 


11 


22 






33 


-5.4090(3) 




91.5(3) 


SJBF3B 






128 


128 


256 


-5.4085(3) 




91.0(3) 


SJBF23B 


4 


8 


128 


128 


268 


-5.4094(3) 




91.8(3) 


PF 












-5.4137(3) 




96.0(3) 


PFBF2B 


11 


22 






33 


-5.4145(3) 




96.8(3) 


PFBF23B 


4 


8 


128 


128 


268 


-5.4152(3) 




97.5(3) 



Est. Exact 



-5.417806 



100.0 
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Table 5.2: Slater- Jastrow (SJ), Pfaffian-Jastrow (PF) and Cl-Jastrow (CI) wave functions 
with backflow (BF) correlations for carbon dimer. The notation is the same as in Table 5.1 1 



Method WF 








N W2 




E[H] 


a 2 [H 2 ] 


E corr [%] 


HF S 


- 


- 


- 


- 


- 


-10.6604 


- 


0.0 


VMC SJ a 


- 


- 


- 


- 


- 


-10.9936(4) 


0.179 


81.7(1) 


SJBF2B 


11 


22 


- 


- 


33 


-11.0012(3) 


0.144 


83.5(1) 


SJBF23B 


4 


8 


128 


128 


268 


-11.0014(2) 


0.141 


83.6(1) 


PF b 


- 


- 


- 


- 


- 


-11.0171(2) 


0.160 


87.4(1) 


PFBF2B 


11 


22 


- 


- 


33 


-11.0223(3) 


0.123 


88.7(1) 


PFBF23B 


4 


8 


128 


128 


268 


-11.0223(2) 


0.128 


88.7(1) 


cr 


- 


- 


- 


- 


- 


-11.0420(4) 


0.112 


93.6(1) 


CIBF2B 


11 


22 


- 


- 


33 


-11.0440(3) 


0.100 


94.0(1) 


DMC SJ 


- 


- 


- 


- 


- 


-11.0227(2) 


- 


88.8(1) 


SJBF2B 


11 


22 


- 


- 


33 


-11.0269(4) 


- 


89.9(1) 


SJBF23B 


4 


8 


128 


128 


268 


-11.0280(3) 


- 


90.1(1) 


PF 


- 


- 


- 


- 


- 


-11.0419(9) 


- 


93.5(2) 


PFBF2B 


11 


22 






33 


-11.0443(6) 




94.1(2) 


PFBF23B 


4 


8 


128 


128 


268 


-11.0447(3) 




94.2(1) 


CI 












-11.0579(5) 




97.5(1) 


CIBF2B 


11 


22 






33 


-11.0580(4) 




97.5(1) 


Est. Exact 












-11.068(5) 




100.0 



a Slater determinant contains PBE DFT orbitals. 

b Same PBE DFT orbitals are used also in PF wave function. 

c Uses natural orbitals with weights of determinants re-optimized in VMC. 



5.4. Conclusions 
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update of all electron positions. In the future, it will be therefore necessary to perform the 
scaling tests to larger systems and compare the overall gain from backflow correlations to 
its overhead cost. 
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Chapter 6 



Summary 

The QMC methodology has proved to be a powerful technique for obtaining the 
ground state properties of fermionic systems. Its only insufficiency comes from the necessity 
to circumvent the fermion-sign problem by the fixed-node approximation. The resulting 
fixed-node errors account for small, but important, fraction of the correlation energies. In 
this dissertation, we have presented developments with direct relevance to elimination or at 
least to improvement of these fixed-node errors. 

In the first part of this dissertation, presented in the Ch. [3j we have analyzed the 
structure and properties of nodes of spin-polarized atomic and molecular fermionic wave 
functions constructed from one-particle orbitals of s, p, d and / symmetries. The study of 
the cases with high symmetries enabled us to find exact nodes for several states with a few 
electrons (p 2 ,7r 2 , and p 3 ). Furthermore, we have used the projection of multi-dimensional 
manifolds into 3D space to study the topologies of the nodes of Hartree-Fock wave functions. 
Finally, we illustrate how the correlation in the accurate CI wave functions, for two specific 
cases of spin-unpolarized states, manifests itself in reducing the nodal structure to only two 
maximal nodal cells. 

In the next chapter, we have proposed a generalized pairing wave function based 
on the Pfaffian functional form. The tests of the Pfaffian pairing wave functions on the set of 
first row atoms and dimers revealed that the Pfaffians were able to capture a large fraction 
of missing correlation energies. We have also explored extensions to linear combinations 
of Pfaffians with good results for atomic systems, but with limited gains in the correlation 
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energies for molecular systems. Further, we have also employed a wave function in the 
form of fully-antisymmetrized independent pairs. We have found that it does not lead to 
additional gains in correlation energy. We conclude that the Pfaffian pairing wave functions 
offer better description of the fermion nodes than the wave functions employing Slater 
determinant (i.e., Hartree-Fock wave functions), and exhibit the correct topology with the 
minimal number of two nodal cells. 

Finally, we have teamed up the Slater, Pfaffian and CI based wave functions with 
inhomogeneous backflow transformation. Our preliminary tests for two chemical systems 
indicate that the backflow correlations reduce the variances of local energy and results 
in some improvement of VMC and DMC energies as well as in better description of the 
fermion nodes. However, these gains come at the price of an additional computational cost, 
justification of which will have to be determined in the near future. 
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Appendix A 



Cutoff-Cusp and Polynomial Pade 



Functions 



The exact electron-electron cusp conditions are satisfied by the choice of following 



cusp-function [see Fig. (A.l)] with variable x = r/r cu t, where r cu t is some cutoff radius and 
7 is the curvature as 

x — x 2 + x 3 /3 1 
1 + 7(3; — x 2 + x 3 /3) 7 + 3 



W*.7) = C ( , - ^ ) (A.1) 



The cusp constant is C = \ for electrons with like and C = \ for electrons with unlike 
spins. 

Polynomial Pade functions for the same variable x = r/r cu t and curvature (5 had 
proved to be excellent choice for describing the electron-electron and electron-nucleus cor- 
relation. In calculations, we use the form 

1 -x 2 (6-8x + 3x 2 ) 

f P oly-Pade{X, P) ~ x + ^ _ ^ + 3^ (^.2) 

The f P oiy-Pade(0) = 1 with derivative f po iy-p a de(®) = ® an d a l so g° es smoothly to zero as 



r cut [see Fig. (A.2)]. These conditions are necessary for preserving cusp conditions 



already fixed by cusp-functions (A.l ) and choice of orbitals or pseudopotentials. 
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Figure A.l: Cutoff-Cusp functions with two different curvatures (7 = 1 and 7 = 10) for like 
(C = |) and unlike (C = |) spins. 
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Figure A. 2: Polynomial Pade functions with curvatures (3 ranging from -0.99 to 100. 
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Appendix B 



Proof of Cayley's Identity 



In order to prove the statement in Eq. (4.7) we will proceed by induction. For 
n = 2 it is true that 



det 



6ia 
-a 12 



pf 



b i2 

-b 12 



pf 



au 
-a 12 



For even n greater than 2, determinant of our matrix of interest can be expanded through 
its cofactors as 








bl2 &13 


&l,n 




-012 


a 2 3 


02,ra 


det 


-Ol3 


-a 23 


«3,n 




_ -«l,n 


—a 2 , n — G3,n • • • 







= ££ 


-ai,kbi,iC{k, 1; 1,/) 





^-a ljfc C(A;,l) 



(B.l) 



The cofactor can be written as 



C{k, 1; 1, Z) = (-l) fc+ ' +1 det [A(k, 1; 1, /)] 



(B.2) 
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where the cofactor matrix is given by 



A(M;1,Z) 



a 23 
-a 23 

-a-2,1 -0,3,1 



-a 2 . 



"03,; 



— 0.k,l 



(B.3) 



At this point we would like to use the induction step and rewrite the determinant cofactor as 



a product of two pfaffians [Cayley's identity Eq. (4.7)]. This would allow us to demonstrate 
that the expansion is identical to the expansion of pfaffians in minors. In order to do so, 
however, we have to shift the k-th column by pair column exchanges, so it becomes the last 
column and, similarly, we have to shift the Z-th row by pair exchanges, so it becomes the 
last row. This involves k pair exchanges of columns and Z pair exchanges or rows and can 
be represented by unitary matrices Uu and U[. It is necessary to invoke these operations so 
that the matrix gets into a form directly amenable for the Cayley's identity, i.e., the matrix 
has to be in a manifestly skew-symmetric form. (The sign change from the row/columns 
exchanges will prove irrelevant as we will show below.) The transformed matrix is given by 



A'(k,l;l,l) = U k A{k,l;l,l)U l 



(B.4) 



and has all zeros on the diagonal with the exception of the last element which is equal to 
—dkj- The last row is given by 



v r =(— a,2,i, • ■ ■ , — a,k-l,i, —ak+i,h ■ ■ ■ 

■ ■ ■ , ~ a i-i,i, «m+i> • • • ' a ',«> ~ a k,l), 
while the last column is given as following 

T i 

v c — \ a 2,k, ■ ■ ■ , Ofc-i,ifc) —a>k,k+x-> ■ ■ ■ 

\T 

■ • • ; — a>k,l-i> —ak,i+i, ■ ■ ■ , —ak,n, — a k,l) ■ 



(B.5) 



(B.6) 



The only non-zero diagonal element —ctki can be eliminated, once we realize that its cofactor 
contains a determinant of a skew-symmetric matrix of odd degree, which always vanishes 
(proof by Jacob [ED]). 
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Now we are ready to perform the induction step, namely to use the property that 
the determinant of a 2(n — 1) x 2(n — 1) matrix can be written as given by the Cayley's 



identity, Eq. (4.7). We obtain 



det[U k A(k, 1; 1, l)U t ] = det[A'(/c, 1; 1,1)] 

= pi[A'(l,k;l,k)}pi[A'(l,l;l,l)}. 



(B.7) 



We can now apply the inverse unitary transformations and shift back the columns (and by 
the skew-symmetry the corresponding rows) in the first pfaffian and, similarly, the rows 
(and corresponding columns) in the second. This enables us to write 

p£[A'(i,k;i,k)] P fL4'(M;U)] 

= pf [U^A(l, k; 1, k)Ui] pf[U k A(l, I; 1, l)U k l ) 

= pf[4(l,fc;l ) fc)]pf[A(l,Z;l J 0] ) (B.8) 



where we have used the identity given by Eq. (4.5d). We can therefore finally write 



C(k, 1; 1, /) = (-l) fe+z+1 pf [A(l, k; 1, fc)]pf L4(l, I; 1, I)] 
= -P c {ai,k)P c {ai,i)y 



(B.9) 



where P c denotes a pfaffian cofactor as defined in (4.4). Therefore, the determinant expan- 



sion in Eq. (B.l ) equals to 

^ -ai,kbi,iC(k, 1; 1, 1) = ^ ai t kh,lP c (ai,k)P c (ai,i) 



k,l 



k.l 

pf[A]pf[B] 



(B.10) 



with matrices A and B defined as in Eq. (4.8 ). This concludes the proof of the more general 



form of the Cayley's identity. Note, if B = A, we trivially obtain well-known formula for 



the square of pfaffian [Eq. (4.5b)]. 
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Appendix C 



Core Pfaffian Algorithms 



C.l Gaussian Elimination with Row-Pivotting Algorithm for 
Pfaffian Value 

int RowPivoting(Array2 <doublevar> & tmp, int i, int n){ 
//row pivoting algorithm used by Pf af f ian_partialpivot 
doublevar big; 
doublevar temp; 
doublevar TINY=le-20; 
Arrayl <doublevar> backup (2*n) ; 
int d=l; 
int k=0; 
big=0.0; 

//find the largest value 
for (int j=i+l; j<2*n; j++){ 
temp=f abs (tmp(i , j)) ; 
if (temp > big){ 
big=temp; 
k=j; 

} 

> 

if (big<TINY){ 

cout «" Singular row in matrix! ! ! "«endl; 
tmp(i,i+l)=TINY; 

> 



C.l. Gaussian Elimination with Row-Pivotting Algorithm for Pfaffian Value 
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if (k!=i+l){ 

//exchange k-th column with 2-nd column; 
for (int j=i; j<2*n; 
backup (j)=tmp(j ,i+l) ; 
tmp(j ,i+l)=tmp(j ,k) ; 
tmp ( j , k) =backup ( j ) ; 

} 

//exchange k-th row with 2-nd row; 
for (int j=i; j<2*n; j++){ 

backup(j)=tmp(i+l, j) ; 

tmp(i+l, j)=tmp(k, j) ; 

tmp (k , j ) =backup ( j ) ; 

} 

d*=-l; //sign change of pfaffian 
return d; 

} 

doublevar Pfaffian_partialpivot (const Array2 <doublevar> & IN){ 
// returns the pfaffian of skew-symmetric matrix IN 
// with partial pivoting 
if (IN.dim[0]°/ 2!=0) return 0.0; 
int n=IN.dim[0]/2; 
Array2 <doublevar> tmp(2*n,2*n) ; 
doublevar PF=1.0; 
doublevar fac; 
for (int i=0;i<2*n;i++) 
for (int l=0;K2*n;l++) 
tmp(i,l)=IN(i,l); 

int d=l; 

for (int i=0;i<2*n;i=i+2){ 

//for given row look for pivoting element 
//exchange if needed 
d*=RowPivoting(tmp, i, n) ; 

for (int j=i+2; j<2*n; j++){ 
fac=-tmp(i, j)/tmp(i,i+l) ; 
for (int k=i+l ;k<2*n;k++){ 

tmp (k , j ) =tmp (k , j ) +f ac*tmp (k , i+ 1 ) ; 

tmp(j ,k)=tmp(j ,k)+f ac*tmp(i+l ,k) ; 

} 



C.2. Algoritm for the Update of Inverse of Pfafnan Matrix 
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} 

PF=PF*tmp(i,i+l) ; 

> 

return PF*d; 

} 



C.2 Algoritm for the Update of Inverse of Pfaffian Matrix 

doublevar UpdatelnversePf af f ianMatrix(Array2 <doublevar> & IN, 

Array 1 <doublevar> & row, 
Array 1 <doublevar> & column, 
int e) 

{ 

//update row and column of skew-symmetric inverse matrix IN 
//the ratio of new/old pfaffians is returned by Column(e) 
int n=in.dim[0]/2; 
for (int i=0;i<2*n;i++){ 

column(i)=0.0; 

for (int j=0; j<2*n; j++) 

column (i)+=row(j ) *IN( j , i) ; 

> 

//to avoid the catastrophe in later division 
if (column(e)==0) 
column (e)=le-20 ; 

//rest is just the new inverse matrix IN 
for (int i=0; i<2*n; i++){ 
if (i==e){ 
IN(i,i)=0.0; 

} 

else { 

IN(e,i)=+IN(e,i)/column(e) ; 
IN(i,e)=-IN(e,i); 

} 

> 

for(int j=0; j<2*n; j++){ 
if (j!=e){ 

for (int k=0;k<2*n;k++){ 
if (k==j) { 
IN(k,k)=0.0; 

} 



C.2. Algoritm for the Update of Inverse of Pfafnan Matrix 



103 



else { 

IN(k, j)-=column(j)*IN(k,e) ; 
IN(j,k)=-IN(k,j) ; 

> 

} 

} 

> 

//return the ratio of pfaffians 
return column (e); 



